Skip to content

Latest commit

 

History

History
257 lines (149 loc) · 15.2 KB

README.md

File metadata and controls

257 lines (149 loc) · 15.2 KB

GramsSky

Table of contents generated with markdown-toc

GramsSky is a simulation of particles coming from spherical (sky) distributions around a detector. Using the options XML file, the user can specify position and energy distributions coming from the inner surface of a "celestial sphere".

Note:

  • If the HepMC3 libraries are not found during the build process, then this package will not be compiled.

  • If the FITSIO and HEALPix libraries are not found during the build process, then any generators involving HEALPix maps won't be compiled into GramsSky.

In the following documentation, names in this format are parameters in the options XML file.

Output

The input and output units for the parameters described in the rest of this documentation are determined by parameters LengthUnit and EnergyUnit. Normally these options are in the <global> section of the options XML file.

The output consists of events in HepMC3 format, which can be used as input to GramsG4.

GramsSky can write files in various formats. The format of the file is assumed to match the file's extension (the part after the final period ".") as follows:

Extension Format
.hepmc3 HepMC3 ASCII
.hepmc2 HepMC2 ASCII
.hpe HEPEVT (ASCII)
.lhef Les Houches Event File
.root HepMC3 ROOT
.treeroot HepMC3 ROOT TTree

If you want to preserve the program options in the output of gramssky as discussed in the Options documentation, you must use the .treeroot format. The program can't save a separate Options ntuple in any other format.

The process of GramsSky generation

The general operation of GramsSky is similar to the General Particle Source available in Geant4. The principle difference is in the random "translation" of the generated particle to simulate that it was generated at "infinity" (or at least several light-years away). Understanding that process is crucial to supplying meaningful options to the program, so it will be addressed first.

The number of events to generate comes from the events parameter in the options file. For each such event:

  • A particle is generated at an original (x,y,z) as determined by the position algorithm selected by the user; see the list of possible algorithms below. Note that the particle's origin is on the surface of a sphere of radius RadiusSphere and center OriginSphere.

    • While nominally OriginSphere would be the center of the detector, the definition of the "center" of the GRAMS detector may depend on the type of study being done: As of 11-Jan-2022, the origin of the z-axis is set to be the bottom of the anode plate (to be partially compatible with LArSoft). In this coordinate system and with the geometry in grams.gdml as of 11-Jan-2022, the center of the LArTPC is (0,0,-14.783); the center of the overall detector is at (0,0,45.217).

  • The particle's kinetic energy E is generated according to the energy algorithm selected by the user; again, see below.

  • The magnitude of the particle's momentum is calculated by p = sqrt(E^2 - m^2), where m is the particle's mass from ROOT's TParticlePDG database based on the value of PrimaryPDG.

  • MapDirection is the direction of "north pole" of the map with respect to the z-axis of the detector coordinate system. This parameter is necessary because the sky maps are normally oriented with (θ,φ)=(90°,0) at the galactic north pole, while a detector will have some other orientation with respect to the earth's surface.

    MapDirection is used to rotate the original (x,y,z) in the celestial sphere coordinate system to (x',y',z') in the detector coordinate system.

  • The momentum direction is taken to be of magnitude p and in the direction of the center of celestial sphere at OriginSphere. Call this (px,py,pz).

  • The direction of (x',y',z') is assumed to be located at a very large distance from the detector. Therefore, all particles from that source appear to be coming from that direction independent of the location within the detector; in other words, the celestial sphere always appears to be at infinity no matter from where you look. To simulate this effect, the source of the particles is taken to be a disc with radius RadiusDisc that is tangent to the celestial sphere at (x',y',z').

    • In the XML file, if RadiusDisc less than or equal to zero, the program will use the value of RadiusSphere. Note that, given the grams.gdml file as of 11-Jan-2022, a radius of 267 cm is sufficient to cover the entire GRAMS outer detector; a radius of 108 cm will cover the LArTPC. (These values come from adding the dimensions of the relevant volume in quadrature to get the maximum length of the volume's diagonal.) Check that the value of OriginSphere is consistent with the RadiusDisc you choose.

      However, if RadiusDisc is larger than those limits, it won't directly affect the simulation. It just means that you're generating more particles whose trajectory won't interact with the detector.

  • The vector with origin (x',y',z') and vector (px,py,pz) is translated along the plane of the disc to a random (r,θ) in the disc's coordinate system. These translated values are what's written as the particle's information in the gramssky output.

    • Note that gramssky does not check to make sure that the translated vertex is within the world volume defined by the GDML file. If a particle's vertex is outside the world volume, gramsg4 will not crash, but it will skip over the particle with a warning message. It's a good idea to verify that the celestial sphere and any potential tangent discs do not lie wholly or partially outside the world volume. One possible fix is to increase the size of the world volume in the GDML file, at a cost of a modest increase in the execution time of gramsg4 as it propagates particles through a larger world.

  • If you are running multiple jobs to generate events, by default they'll all run with the same random number seed; i.e., in the options XML file there is a parameter rngseed which is set to -1 by default. To generate a different set of events for each job, you will want to vary the seed for each job.

    For example, if the job has a unique process ID in the variable ${Process}, then you probably want something like this:

    ./gramssky --rngseed ${Process}

Position generators

These are the available generators for (x,y,z) as 11-Jan-2022. If option PositionGeneration has the value:

"Point"

The user must supply the option PointSource with a vector of values (x,y,z). All the particles are generated with an initial vertex at that point, then undergo the geometric transformations described above.

  • The position (x,y,z) will be translated to the surface of the celestial sphere automatically. There's no need to do elaborate calculations to make sure (x,y,z) is exactly RadiusSphere from OriginSphere; the program takes care of that.

"Iso"

The particle position is generated isotropically on the inner surface of the celestial sphere. The controlling options are:

  • ThetaMinMax: the value is a two-element vector with minimum and maximum θ.
  • PhiMinMax: the value is a two-element vector with minimum and maximum φ.

Units are radians. Setting ThetaMinMax to (0,1.571) and PhiMinMax to (0,6.283) will generate particle isotropically in a hemisphere above the detector (subject to the geometry transforms described above).

Energy generators

These are the available generators for E as 11-Jan-2022. If option EnergyGeneration has the value:

"Fixed"

The value of E comes from the option FixedEnergy.

"Flat"

E is generated uniformly between the values of parameters EnergyMin and EnergyMax.

"Gaus"

The gaussian distribution is of the form:

where μ is the mean of the distribution and σ is the width. E is generated with μ given by GausMean and σ given by GausWidth.

  • Note that the limits in parameters EnergyMin and EnergyMax still apply to this generator. This is to keep the value of E from going negative, which would cause problems in both gramssky and gramsg4.

"BlackBody"

The black-body radiation function is of the form:

where kT is the "radiation temperature". E is generated according to a black-body distribution with kT given by parameter RadTemp.

  • Again, the limits in parameters EnergyMin and EnergyMax still apply to this generator.
  • The units of kT must be the same as that of EnergyMin and EnergyMax; i.e., the value of the global parameter EnergyUnit.

"PowerLaw"

The power-law function is of the form:

where

  • N is a normalization
  • Eref is the "reference energy"
  • α is the "photon index"

E is generated according to a power-law distribution with α given by parameter PhotonIndex.

  • The limits in parameters EnergyMin and EnergyMax still apply to this generator.

"Hist"

The program expects two parameters: HistFile with the name of a ROOT file, and HistName with the name of a histogram within that file.

  • Both HistFile and HistName can contain path specifiers for directories either on the computer system or within the ROOT file.

  • For this generator, the values of EnergyMin and EnergyMax are ignored. Instead, the energy limits effectively come from the bin limits of the histogram.

Combined position and energy generators

As of Feb-2022, all the generators in this category make use of the FITSIO and HEALPix libraries.

  • FITS is a file format intended for both images and multi-dimensional data.

  • HEALPix is a pixelization for evenly subdividing a sphere.
Courtesy NASA/JPL-Caltech

  • Credit to Naomi Tsuji and Hiroki Yoneda, who provided me with the code and files to incorporate the following into GramsSky.

  • As noted above, if the healpix_cxx libraries are not installed on your system, the following generators will not be compiled into GramsSky.

For the following options, the value of EnergyGeneration is ignored. If option PositionGeneration has the value:

"MapPowerLaw"

This method uses three HEALPix maps, one for each parameter in a power-law distribution:

The procedure is to randomly select a pixel, then randomly generate the energy according to the power-law distribution at that position.

This approach is intended as a simple simulation for stellar sources.

The parameters for "MapPowerLaw" are:

  • "MapPowerLawFile" = the name of the HEALPix file containing the power-law maps.

  • "MapPowerLawHDU" = the HDU for the HEALPix maps within the file.

  • "MapPowerLawColumnNorm" = the "column number" of the map for parameter N within the HDU.

  • "MapPowerLawColumnIndex" = the "column number" of the map for parameter α within the HDU.

  • "MapPowerLawColumnEref" = the "column number" of the map for parameter Eref within the HDU.

  • The limits in parameters EnergyMin and EnergyMax are applied to the power-law distribution within each pixel.

Note that this HDU/column structure may not be permanent, depending on the evolution of the process as determined by Naomi Tsuji and Hiroki Yoneda.

"MapEnergyBands"

This method uses a series of HEALPix maps, one for each of an increasing set of energies.

The procedure is to randomly select an energy band, then randomly select a pixel from the flux in that band, and finally generate an energy distribution for that particular pixel from a power-law distribution.

This approach is intended to simulate the diffuse sky background.

The parameters for "MapEnergyBands" are:

  • "MapEnergyBandsFile" = the name of the HEALPix file containing the energy-band maps.

  • "MapEnergyBandsHDU" = the HDU for the HEALPix maps within the file.

  • "MapNumberEnergyBandsKey" = within the HDU, this is a key whose value is the number of energy-band maps in the file.

  • "MapEnergyBandsPrefix" = The key for each map is formed by this string, suffixed by a number. For example, if MapNumberEnergyBandsKey is "NMAP" and MapEnergyBandsPrefix is "ENE", the individual maps have keys "ENEnn" where nn is 1 through NMAP.

  • Only those energy bands with energies beween EnergyMin and EnergyMax are used in the above procedure.

Note that this HDU/column structure may not be permanent, depending on the evolution of the process as determined by Naomi Tsuji and Hiroki Yoneda.