This is a Geant4 [1] [2] [3] simulation of the Upstream Target Room (UTR) at the High Intensity γ-ray Source (HIγS) facility [4], located at the Duke Free Electron Laser Laboratory (DFELL) of the Triangle University Nuclear Laboratory (TUNL) at Duke University, Durham, NC, USA.
-
2.1 Geometry
2.3 Event Generation
2.4 Physics
-
3.1 Dependencies
3.2 Compilation
This section is an overview or checklist for a user who wants to implement his own geometry and do simulations.
See section 3 Installation
See section 4 Usage and Visualization
- Create a directory in
DetectorConstruction. - Implement a geometry in a
DetectorConstruction.hhandDetectorConstruction.ccfile in this directory. - Configure
utrto use the new geometry using the cmake build option (see sections 2.1 Geometry or 3 Installation)
- In
DetectorConstruction.cc, set volumes as sensitive detectors. There is the choice between 3 different detector types that record different information about particles (see section 2.2 Sensitive Detectors). - In
src/ActionInitialization.cc, set the quantities that should be written to the output file (see section 2.2 Sensitive Detectors).
Configure utr to use either the G4GeneralParticleSource or the AngularDistributionGenerator using the CMake build option (see sections 2.3 Event Generation or 3 Installation)
In src/Physics.cc activate or deactivate processes, choose between different models, or include certain particles (see section 2.4 Physics)
In utr.cc, set the random number seed explicitly to get deterministic results. It is not clear whether this works in multithreaded mode. See section 2.5 Random Number Engine
Create a macro file that contains instructions for the primary generator (see section 2.3 Event Generation).
Do not simulate more than 2^32 ~ 2 billion particles using /run/beamOn, since this causes an overflow in the random number seed, giving you in principle the same results over and over again. Execute the same simulation multiple times instead.
Consider creating a macro like loop.mac to loop over variables in your macro file.
(See section 4 Usage and Visualization)
See section 5 Output Processing
The DetectorConstruction.cc has to be adapted by the user for their experiment. Several pre-defined geometries from various experiments already exist in the DetectorConstruction/ directory. The most recent geometries are from the experiments on 82Kr/82Se and 120Sn/82Kr. Those have implementations of the most relevant parts of the collimator room and the utr (γ³ setup [5] and second setup).
By default, utr will use the geometry in DetectorConstruction/82Se_82Kr_726_743 for no particular reason. To use the geometry in "DetectorConstruction/DetectorConstructionXY", set the corresponding CMake build variable when building the source code (see also 3 Installation).
$ cmake -DDETECTOR_CONSTRUCTION="DetectorConstructionXY" .Several pre-defined classes exist to simplify the geometry construction which are explained in the following.
Classes for several different detectors exist. In all of those, a G4LogicalVolume that contains all the parts of a detector is implemented which is returned by the class' Get_Logical() method. Furthermore, each detector class can return its mother volume's length and radius.
To place a detector DetectorXY in the geometry, create an instance of the detector class in DetectorConstruction.cc, get the logical volume and place it in the geometry (using Get_Length() and Get_Radius() if necessary):
DetectorXY* detectorXY_Instance = new DetectorXY("DetectorXY_Name");
G4LogicalVolume* DetectorXY_Logical = detectorXY_Instance->Get_Logical();
new G4PVPlacement( ... )
The name "DetectorXY_Name" is the name of the logical volume of the detector crystal which can be used to register it as a sensitive detector.
At the moment, the following detectors are implemented:
- Cologne LaBr, Saint Gobain BrilLanCe 380 1.5x1.5"
- Darmstadt LaBr, Saint Gobain BrilLanCe 380, 3x3"
- Duke 60% HPGe "1" (Ortec serial number 36-TN31061A)
- Duke 60% HPGe "2" (Ortec serial number 36-TN30986A)
- Duke 60% HPGe "3" (Ortec serial number 36-TN40663A)
- Duke 60% HPGe "4" (Ortec serial number 36-TN21033A)
- Duke 55% HPGe (Ortec serial number 4-TN21638A)
- Duke 55% HPGe (Ortec serial number 4-TN31524A)
- Duke 120% HPGe "Zero degree detector" (Ortec serial number 33-P40383A)
- Darmstadt HPGe "1" (Canberra serial number 10PC473156-A)
- Darmstadt HPGe "2" (Eurisys Mesures serial number 10PC447589-A)
- Darmstadt Clover "Polarimeter" (Eurisys Mesures serial number 10PC447590-a)
For maximum flexibility, the shielding of the setup can be constructed brick by brick. To avoid the G4Solid->G4LogicalVolume->G4PhysicalVolume procedure each time one would like to place a standardized brick, a lot of them are predefined as classes in Bricks.hh.
Once instantiated in DetectorConstruction.cc, bricks can be placed inside the G4LogicalVolume which was defined to be their mother volume via their constructor. To place a brick, use the Put(x, y, z, angle_x, angle_y, angle_z) method in which the coordinates and rotation angles around the coordinate axes can be specified.
Bricks are assumed to be cuboid objects, i.e. they can have 3 different side lengths. In Bricks.hh, the convention is that the long side points in z-direction, the medium side in x-direction and the short side in y-direction, if they can be distinguished. The respective lengths can be accessed via the member variables L, M and S.
Similar to bricks, filters and filter cases in front of detectors are implemented in Filters.hh and can be placed using their Put() methods.
The intent of this was to give the user an overview which filters are really there at the UTR. Sometimes, the documentation of experiments only mentions "thin", "medium" and "thick" filters, but no actual dimensions. Providing fixed filter types hopefully helps with such issues.
Complicated targets can be implemented in Targets.hh. The placement in DetectorConstruction.cc works analog to the placement of detectors. Relevant properties of the targets can be made accessible by implementing Get() methods.
Information about the simulated particles is recorded by instances of the G4VSensitiveDetector class. Any unique logical volume can be declared a sensitive detector.
Unique means volumes with a unique logical volume name. This precaution is here because all bricks and filters of the same type from the previous sections have the same logical volume names. Making one of those a sensitive detector might yield unexpected results.
Any time a particle executes a step inside a G4VSensitiveDetector object, its ProcessHits routine will access information of the step. This way, live information about a particle can be accessed. Note that a "hit" in the GEANT4 sense does not necessarily imply an interaction with the sensitive detector. Any volume crossing is also a hit. Therefore, also non-interacting geantinos can generate hits, making them a nice tool to explore to explore the geometry, measure solid-angle coverage, etc. ... After a complete event, a collection of all hits inside a given volume will be accessible via its HitsCollection. This way, cumulative information like the energy deposition inside the volume can be accessed.
Three types of sensitive detectors are implemented at the moment:
- EnergyDepositionSD Records the total energy deposition by any particle per single event inside the sensitive detector.
- ParticleSD Records the first step of any particle inside the sensitive detector.
- SecondarySD Records the first step of any secondary particle inside the sensitive detector.
No matter which type of sensitive detector is chosen, the simulation output will be a ROOT tree with a user-defined subset (see section 2.6 Output File Format) of the following 10 branches:
- ekin Kinetic energy (in MeV) of the particle at the time of its first hit of the sensitive detector.
- edep Total energy deposition (in MeV) of the event (EnergyDepositionSD) OR energy deposition of the first hit of the sensitive detector (ParticleSD, SecondarySD)
- particle Type of the particle whose first hit of the sensitive detector it was (ParticleSD, SecondarySD) OR type of the very first particle in this event that hit the sensitive detector (EnergyDepositionSD). The type of the particle is encoded in the Monte Carlo Particle Numbering Scheme.
- volume User-defined identifier of the sensitive detector that recorded this set of data.
- x/y/z Coordinates (in mm) of the first hit of the sensitive detector by a particle (ParticleSD, SecondarySD) OR coordinates of the first hit by the first particle in this event that hit the sensitive detector (EnergyDepositionSD)
- vx/vy/vz Momentum (in MeV/c) of the particle at the position of the first hit of the sensitive detector (ParticleSD, SecondarySD) OR momentum of the first particle hitting the sensitive detector in this event at the position of its first hit (EnergyDepositionSD).
The meaning of the columns sometimes changes with the choice of the sensitive detector.
To make a volume whose logical volume name is "Logic_Name" a sensitive detector of type XYSD, add the following lines in DetectorConstruction::ConstructSDandField(){}
XYSD *xySD = new XYSD("SD_Name", "SD_HitsCollection_Name");
G4SDManager::GetSDMpointer()->AddNewDetector(xySD);
xySD->SetDetectorID(volume); // volume is the detector ID which will be displayed in the ROOT output file
SetSensitiveDetector("Logic_Name", xySD, true);
The following examples illustrates how the different sensitive detectors work.
In the figure, a photon (orange, sinusoidal line) and several electrons (blue arrow) are shown which are part of a single event. Each hit in the sensitive detector (black circle) has a label that contains the particle number and the number of the hit. The history of the event is as follows:
The primary particle 1 (a photon with an energy of 2.5 MeV) enters the sensitive detector at (0,3). It is Compton-scattered at (2,5) and transfers 1.0 MeV to an electron (particle 2). The electron slowly loses its kinetic energy again in scattering processes at (4,8), (6,9) and (7,8) inside the sensitive detector which do not create new secondary particles. In the meantime, particle 1 travels to (6,1) and creates an e-/e+ pair with its remaining 1.5 MeV. Each lepton (3 and 4) gets an initial energy of (1.500 - 1.022)/2 = 0.239 MeV. Lepton 3 loses its kinetic energy inside the sensitive detector, Lepton 4 leaves the sensitive detector.
The following information would be recorded if the volume was a ...
- ... ParticleSD
| ekin | edep | particle | volume | x | y | vx | vy |
|---|---|---|---|---|---|---|---|
| 2.5 | 0. | 22 | 0 | 0 | 3 | 2 | 2 |
| 1.0 | 0. | 11 | 0 | 2 | 5 | 2 | 3 |
| 0.239 | 0. | 11 | 0 | 6 | 1 | 3 | 3 |
| 0.239 | 0. | 11 | 0 | 6 | 1 | -3 | -3 |
- ... SecondarySD
| ekin | edep | particle | volume | x | y | vx | vy |
|---|---|---|---|---|---|---|---|
| 1.0 | 0. | 11 | 0 | 2 | 5 | 2 | 3 |
| 0.239 | 0. | 11 | 0 | 6 | 1 | 3 | 3 |
| 0.239 | 0. | 11 | 0 | 6 | 1 | -3 | -3 |
- ... EnergyDepositionSD
| ekin | edep | particle | volume | x | y | vx | vy |
|---|---|---|---|---|---|---|---|
| 2.5 | 1.989 | 22 | 0 | 0 | 3 | 2 | 2 |
For simplicity, the momentum vector is given as the vector of the depicted arrow/wave. The real values would be
vx_real = 1/sqrt(vx^2 + vy^2)*cos(arctan(y/x))*ekin/c
vy_real = 1/sqrt(vx^2 + vy^2)*sin(arctan(y/x))*ekin/c
Event generation is done by classes derived from the G4VUserPrimaryGeneratorAction. In the following, the two existing event generators are described.
By default, utr uses the Geant4 standard G4GeneralParticleSource. To use the AngularDistributionGenerator of utr, which implements angular distributions for Nuclear Resonance Fluorescence (NRF) applications, unset the USE_GPS option when building the source code (see also 3 Installation):
$ cmake -DUSE_GPS=OFF .This event generator uses the G4GeneralParticleSource (GPS) whose parameters can be controlled by a macro file. For further information, see the GPS documentation.
The AngularDistributionGenerator generates monoenergetic particles that originate in a set of G4PhysicalVolumes of the DetectorConstruction and have a certain angular distribution.
The physical volumes (the "sources") and angular distribution can have arbitrary shapes. Starting positions and momentum directions are created using a Monte-Carlo method: Given a(n)
- Cuboid in 3 dimensions ("container volume") that completely contains the source volume
- Angular distribution W(θ, φ)
AngularDistributionGenerator generates uniform random
- Positions
(random_x, random_y, random_z) - Tuples
(random_θ, random_φ, random_W)
until
- The position
(random_x, random_y, random_z)is inside one of the source physical volumes random_W <= W(random_θ, random_φ)
If both conditions are fulfilled, a particle is emitted from (random_x, random_y, random_z) in the direction (θ, φ).
The process of finding a starting point is depicted in 2 dimensions in the figure below. The yellow source volume consists of two disconnected parts. It is placed inside a dashed box which is the container volume. Random points inside the container volume are either black or red, depending on whether they are inside or outside the source volumes. Only from the black points, particles are emitted (in random directions).
The process of finding a starting vector is shown in one dimension (W is only dependent on θ) in in the figure below. First, a random value random_θ for θ with a uniform random distribution on a sphere. Note that this is not the same as a uniform distribution of values between 0 and π for θ. Then, a unform random number between 0 and an upper limit W_max is drawn. If the value W_max is lower than W(random_θ) (black points), then a particle will be emitted at that angle.
Clearly, the algorithm works well if
- the cuboid approximates the shape of the sources well and wraps it tightly
- the angular distribution varies smoothly in the
(θ, φ)plane and W_max is very close to the actual maximum value ofW
The maximum number of randomly sampled points is hard-coded AngularDistributionGenerator.cc where it says:
MAX_TRIES_POSITION = 10000
MAX_TRIES_MOMENTUM = 10000
AngularDistributionGenerator can do a self-check before the actual simulation in which it creates MAX_TRIES_XY points and evaluates how many of the were valid. From this, the probability p of not hitting one of the source volumes / angular distribution can be estimated. An individual check is done for each source volume. If p * N >~ 1, where N is the number of particles to be simulated, the algorithm will very probably fail once in a while so try increasing MAX_TRIES_XY. A typical output of the self-check for the position generator looks like:
G4WT0 > ========================================================================
G4WT0 > Checking Monte-Carlo position generator with 100 3D points ...
G4WT0 > Check finished. Of 100 random 3D points, 82 were inside volume 'Se82_Target' ( 82 % )
G4WT0 > Probability of failure: pow( 0.18, 100 ) = 3.36706e-73 %
G4WT0 > ========================================================================
The momentum generator contains one more piece of information. This is because the value of W_max that was introduced above has been arbitrarily set to 3 in AngularDistributionGenerator.hh, which should be a sensible choice for most angular distributions. However, it may be that W(θ, φ) > 3 for some distribution. This can also be detected by the self-check with a Monte-Carlo method. For each of the MAX_TRIES_MOMENTUM tries, utr will also check whether the inequality W_max <= W(random_θ, random_φ) holds. If yes, this could be a hint that the value of W_max is too low and should be increased manually.
The corresponding message would look like
G4WT0 > In 5624 out of 10000 cases (56.24 % )
W(random_theta, random_phi) == MAX_W == 1 was valid.
This may mean that MAX_W is set too low and
the angular distribution is truncated.
If everything is okay, it will display
G4WT0 > MAX_W == 3 seems to be high enough.
However, using W_max = 3 might still be inefficient, so it is good to know the minimum and maximum value of the angular distribution in advance. For debugging of the angular distribution, a unit test exists in the repository, which is especially recommended if one defines new angular distributions.
To change parameters of the AngularDistributionGenerator, an AngularDistributionMessenger has been implemented that makes the following macro commands available:
/ang/particle NAMESet the particle type by name. See the GEANT4 documentation for a list of possible particles./ang/energy ENERGY UNITSet the kinetic energy of the particle/ang/nstates NSTATESDefine the number of states for the angular distribution (see below)/ang/stateN VALUEDefine state number N (where N is in {1, 2, 3, 4})/ang/deltaN1N2 VALUEDefine the multipole mixing ratio for the transition between states N1 and N2/ang/sourceX VALUE UNITAlong withsourceYandsourceZ, defines the position of the container box/ang/sourceDX VALUE UNITAlong withsourceDYandsourceDZ, defines the dimensions of the container box./ang/sourcePV VALUEEnter the name of a physical volume that should act as a source. To add more physical volumes, call/ang/sourcePVmultiple times with different arguments (about using multiple sources, see also the caveat at the end of this section).
The container volume's inside will be the interval [X - DX/2, X + DX/2], [Y - DY/2, Y + DY/2] and [Z - DZ/2, Z + DZ/2].
The identifiers of the angular distribution are given to the simulation as an array of numbers states = {state1, state2, ...} whose length NSTATES can to be specified by the user.
For "real" NRF angular distributions, this array of numbers will be the spins of the excited states in a cascade, with the parity indicated by the sign of the numbers.
Due to parity symmetry, an inversion of all the parities in a cascade will result in the same angular distribution. For example, the transition 0+ -> 1- -> 0+ has the same angular distribution as 0- -> 1+ -> 0-. This is indicated by the notation "+-" and "-+" in the list below (the example is listed as "0+- -> 1-+ -> 0+-").
It was chosen to represent the parity of a state as the sign of the spin quantum number. Unfortunately, this makes it impossible to represent 0- states, because the number "-0" is the same as "+0". Therefore, the value "-0.1" represents a 0- state.
At the moment, the following distributions are implemented:
states == {0.1, 0.1, 0.1}Wildcard for test distributionsstates == {0., 0., 0.}Isotropic distributionstates == {0., 1., 0.} or {-0.1, -1., -0.1}0+- -> 1+- -> 0+-states == {0., -1., 0.} or {-0.1, 1., -0.1}0+- -> 1-+ -> 0+-states == {0., 2., 0.} or {-0.1, -2., -0.1}0+- -> 2+- -> 0+-states == {0., 2., 2.} or {-0.1, -2., -2.}0+- -> 2+- -> 2+-states == {0., 1., 2.} or {-0.1, -1., -2.}0+- -> 1+- -> 2+-states == {0., -1., 2.} or {-0.1, 1., -2.}0+- -> 1-+ -> 2+-states == {1.5, -2.5, 1.5} or {-1.5, 2.5, -1.5}3/2+- -> 5/2-+ -> 3/2+-states == {1.5, 2.5, 1.5} or {-1.5, -2.5, -1.5}3/2+- -> 5/2+- -> 3/2+-states == {1.5, 1.5, 1.5} or {-1.5, -1.5, -1.5}3/2+- -> 3/2+- -> 3/2+-states == {1.5, -1.5, 1.5} or {-1.5, 1.5, -1.5}3/2+- -> 3/2-+ -> 3/2+-states == {0.5, 1.5, 0.5} or {-0.5, -1.5, -0.5}1/2+- -> 3/2+- -> 1/2+-states == {0.5, -1.5, 0.5} or {-0.5, 1.5, -0.5}1/2+- -> 3/2-+ -> 1/2+-states == {2.5, 1.5, 2.5} or {-2.5, -1.5, -2.5}5/2+- -> 3/2+- -> 5/2+-states == {2.5, -1.5, 2.5} or {-2.5, 1.5, -2.5}5/2+- -> 3/2-+ -> 5/2+-
Finding the correct dimensions of the container box might need visualization. Try placing a G4Box with the desired dimensions at the desired position in the geometry and see whether it encloses the source volume completely and as close as possible. In fact, most of the predefined geometries have such a G4Box called AuxBox which is commented out. After using it to determine the size and position of the container box, remember to comment out the code again.
When using multiple sources, be aware that AngularDistributionGenerator samples the points with a uniform random distribution inside the container volume. How many particles are emitted from a certain part of the source will only depend on its volume.
This is not always desirable. Imagine the following example: The goal is to simulate a beam on two disconnected targets. The first target has twice the volume of the second target, but the second target has a four times larger density. That means, the reaction would occur approximately twice as often in the second target than in the first. Implementing the targets like this in AngularDistributionGenerator would not give realistic results because of the weighting by volume. In a case like this, one would rather simulate the individual parts of the target and compute a weighted sum of the results.
The physics processes are defined in /src/Physics.cc. In the Physics::ConstructProcess() method, physics lists for groups of particles can be activated by uncommenting the desired physics list. If you are working on a slow machine and you just want to visualize the geometry, it can be advantageous to switch off all physics.
At the moment, physics processes for gammas, electrons/positrons, neutrons and other charged particles are available. Check Physics.cc to see which processes are defined. A description of the processes can be found in the GEANT4 Physics Reference Manual
The physics list is optimized for NRF applications and largely inspired by the GEANT4 examples TestEm7 and Hadr04.
Mandatory for all particles. Without this, the simulation will do almost nothing.
Uses the Penelope OR Livermore low-energy electromagnetic physics framework for gammas, electrons and positrons. Only one of them can be active at the same time.
Both frameworks aim at a more precise description of these processes than the standard G4 physics which are optimized for particle physics applications at CERN.
A specialized neutron physics framework for low energies.
Standard processes for any other charged particle that might appear in the simulation. The probability for this in NRF applications is low, so this list is not very detailed.
In utr.cc, the random number engine's seed is set by using the current CPU time, making it a "real" random generator.
G4Random::setTheEngine(new CLHEP::RanecuEngine);
// 'Real' random results
time_t timer;
G4Random::setTheSeed(time(&timer));
If you want deterministic results for some reason, comment these lines out and uncomment
G4Random::setTheSeed(long n);
Re-compile afterwards. After this change, every restart of the simulation with unchanged code will yield the same results.
In section 2.2 Sensitive Detectors the format of the ROOT output file was already introduced. The possible branches are
- ekin
- edep
- particle
- volume
- x/y/z
- vx/vy/vz
The user can specify in ActionInitialization.cc which of these quantities should be written to the ROOT file, to avoid creating unnecessarily large files. When the RunAction is initialized in ActionInitialization.cc, set the value of the corresponding output flag to 0 if this quantitiy should not be written to the ROOT file.
For example, the code
output_flags[EKIN] = 0;
output_flags[EDEP] = 1;
output_flags[PARTICLE] = 1;
output_flags[VOLUME] = 1;
output_flags[POSX] = 0;
output_flags[POSY] = 0;
output_flags[POSZ] = 0;
output_flags[MOMX] = 0;
output_flags[MOMY] = 0;
output_flags[MOMZ] = 0;
would only write the energy deposition, the particle type and the senstive detector volume.
These instructions will get you a copy of the simulation running.
To build and run the simulation, the following dependencies are required:
- Geant4 (tested with version 10.04)
- CMake (build)
- Make (build)
Furthermore, to use the analysis scripts:
- ROOT (tested with version 5.34/36)
Optional components:
- Qt as a visualization driver. Compile Geant4 with the
GEANT4_USE_QToption (tested with Qt4) - ccmake UI for CMake which gives an overview of available build options
The most simple way to build and compile the simulation is
$ cmake .
$ makeThis will compile the simulation using a default geometry (see 2.1 Geometry) and the GeneralParticleSource as a primary generator (see 2.3 Event Generation).
To use another geometry in the directory DetectorConstruction/DetectorConstructionXY, set the DETECTOR_CONSTRUCTION build variable when executing the cmake command:
$ cmake -DDETECTOR_CONSTRUCTION="DetectorConstructionXY" .To use the AngularDistributionGenerator instead of the GeneralParticleSource, set the USE_GPS variable
$ cmake -DUSE_GPS=OFF .Any of the commands above will create the utr binary in the top directory.
In case CMake complains that it cannot find Geant4Config.cmake or geant4-config.cmake, set the CMAKE_INSTALL_PREFIX variable to the directory where GEANT4 has been installed:
$ cmake -DCMAKE_INSTALL_PREFIX=/PATH/TO/G4/INSTALL .
$ makeThe compiled utr binary can be run with different arguments. To get an overview, type
$ ./utr --helpImportant optional arguments besides --help are:
$ ./utr -m MACROFILEExecutes macro file MACROFILE
$ ./utr -t NTHREADSSets the number of threads in multithreaded mode (default: 1)
$ ./utr -o OUTPUTDIRSets the output directory of utr where the ROOT files will be placed.
Running utr without any argument will launch a UI session where macro commands can be entered. It should also automatically execute the macro file init_vis.mac, which visualizes the geometry.
If this does not work, or to execute any other macro file MACROFILE, type
/control/execute MACROFILEin the UI session. In this mode, it is important to know how the angles are defined, to be able to set the view using the /vis/viewer/set/viewpointThetaPhi command:
/vis/viewer/set/viewpointThetaPhi 180 0 degviews the setup in beam direction,
/vis/viewer/set/viewpointThetaPhi 0 0 degviews it against the beam direction. To view the geometry from above and below, use
/vis/viewer/set/viewpointThetaPhi 270 270 degand
/vis/viewer/set/viewpointThetaPhi 270 90 degrespectively.
It is also possible to create 3D visualization files that can be viewed by an external viewer like Blender (the title picture was made in Blender, for example). The macro vrml.mac shows how to create a .wrl file. Run it in UI mode with
/control/execute vrml.macThe directory OutputProcessing contains some sample ROOT scripts that can be adapted by the user to process their simulation output. Executing
$ cd OutputProcessing/
$ makein this directory should compile all the scripts and move the executables to the utr/ directory.
The compilation may fail if the ROOTSYS environment variable is not set on your system.
Executables can be run like
$ ./EXECUTABLENAME {ARGUMENTS}The executables can be removed using the command
$ make cleanin OutputProcessing
RootToTxt converts a ROOT output file (TFile) containing an n-tuple of data (a TTree with TBranch objects) to a simple text file with the same content. If you want to convert a ROOT file ROOTFILE, type
$ ./rootToTxt ROOTFILEThe ROOT file can have a TTree with an arbitrary name and an arbitrary number of TBranch objects. The output text file has the same name as the ROOT file but with a ".txt" suffix. Be aware that conversion into text files increases the file size.
getHistogram sorts the data from multiple output files (for example, several threads of the same simulation) into a ROOT histogram and saves the histogram to a new file.
Executing
$ ./getHistogram --help
Usage: getHistogram [OPTION...] Create histograms from a list of events
GetHistogram
-b BIN Number of energy bin whose value should be
displayed
-m MULTIPLICITY Particle multiplicity
-o OUTPUTFILENAME Output file name
-p PATTERN1 File name pattern 1
-q PATTERN2 File name pattern 2
-t TREENAME Name of tree
-v Verbose mode (does not silence -b option)
-?, --help Give this help list
--usage Give a short usage messageshows how to use the script. The meaning of the arguments to the options is:
- TREE: Name of the ROOT tree (Default: TREE=="utr") in all of the output files.
- PATTERN1 and 2: Two strings that identify the set of files to be merged. See also the example below. (Default: PATTERN1=="utr", PATTERN2==".root")
- OUTPUTFILE: Name of the output file that contains the histograms. (Default: OUTPUTFILENAME=="hist.root")
- MULTIPLICITY: Determines how many events should be accumulated before adding information to the histogram. This can be used, for example, to simulate higher multiplicity events in a detector: Imagine two photons with energies of 511 keV hit a detector and deposit all their energy. However, the two events cannot be distinguished by the detector due to pileup, so a single event with an energy of 1022 keV will be added to the spectrum. Similarly, Geant4 simulates events by event. In order to simulate pileup of n events, set MULTIPLICITY==n. (Default: MULTIPLICITY==1)
- BIN: Number of the histogram bin that should be printed to the screen while executing
getHistogram. This option was introduced because often, one is only interested in the content of a special bin in the histograms (for example the full-energy peak). If the histograms are defined such that bin3000contains the events with an energy deposition between2.9995 MeVand3.0005 MeVand so on, so there is an easy correspondence between bin number and energy. (Default: BIN==-1)
A short example: The typical output of two different simulations on 2 threads each are the files
$ ls
master.root
utr0_t0.root
utr0_t1.root
utr1_t0.root
utr1_t1.root
(The master.root file may be there but it contains no data. It is just an artifact of the multithreaded mode.)
Assuming you would like to merge both threads of utr0 and write the accumulated histogram to hist.root, you execute
$ ./getHistogram -t utr -p utr0 -q .root -o hist.rootThis will create an output file called hist.root with the desired histogram.
If, instead, you wanted to merge all ROOT files, do
$ ./getHistogram -t utr -p utr -q .root -o hist.rootIn this simple example, already the first pattern uniquely identifies the files you want to merge. However, imagine there was a file utr0.txt in the same directory, then the second pattern could be used to exclude this.
Yet another possibility is that you would like to merge both threads of utr0 into one histogram, and both threads of utr1 into another one. The shell script loopGetHistogram.sh can be used for this task. Executing
$ ./loopGetHistogram 0 1 utr utrwould do just what is described above, creating the output files utr0.root and utr1.root.
The user has to hard-code the desired histograms in OutputProcessing/GetHistogram.cpp, i.e. edit the lines from
//
// START OF USER-DEFINED OUTPUT
//and
//
// END OF USER-DEFINED OUTPUT
//The standard getHistogram script that is included in this repository creates histograms (0 to 10 MeV, 10^4 bins) for the energy deposition by photons in the volumes with numbers 1 to 9 (about setting the number of a sensitive volume, see 2.2 Sensitive Detectors). Setting a higher multiplicity means that the energy deposition of several simulated events is summed up.
A direct follow-up to getHistogram, HistogramToTxt.cpp takes a ROOT file that contains only 1D histograms (TH1 objects) and converts each histogram to a single text file. The script is used as follows:
$ ./histogramToTxt FILENAME.rootWhere FILENAME is the name of the ROOT file that contains TH1 objects. histogramToTxt outputs text files with the naming pattern
HISTNAME_FILENAME.txtwhere HISTNAME is the name of the TH1 object and FILENAME the same as above.
The shell script loopHistogramToTxt.sh shows how to loop the script over a large number of files.
For testing the AngularDistributionGenerator, a dedicated geometry can be found in /DetectorConstruction/AngularDistributionGenerator_Test/ and a macro file and output processing script are located in /unit_tests/AngularDistributionGenerator/. The test geometry consists of a very small spherical particle source surrounded by a large hollow sphere which acts as a ParticleSD. Geantino particles emitted by this source and detected by the hollow sphere should have the desired angular distribution. This test was originally implemented to test the built-in angular distributions which are manually coded from the output of a computer algebra program, and to get a feeling of how large the value of W_max has to be.
In order to use the unit test, the following things have to be prepared:
- Build the correct geometry by setting the build variable
DETECTOR_CONSTRUCTIONtoAngularDistributionGenerator_Test(see 3.2 Compilation). - Ensure that the momentum vectors of the particles are written to a ROOT file by enabling the correct flags in
ActionInitialization.hh(see 2.6 Output File Format) - Run a simulation with geantinos and the desired angular distribution. A sample macro is
/unit_tests/AngularDistributionGenerator/angdist_test.mac - Compile the analysis script
AngularDistributionGenerator_Test.cppin the same directory by typingmake. It can be executed with a number of simulation output files similar toGetHistogram(see 5. Output Processing).
When executed, the analysis script creates a file hist.root that contains a histogram in the (θ, φ) plane which shows the frequency of particles with a certain direction. This can already be visually compared to the shape of the continuous function that was used as an input. An example for a 0^+ -> 1^+ -> 0^+ cascade is shown below:
With the bad choice W_max == 1, the same analysis would result in the following histogram which is cleary different:
In order to have a more quantitative analysis, the histogram can also be fit with the input distribution as it is implemented in src/AngularDistribution.cc. For this purpose, AngularDistribution.hh is already included in AngularDistributionGenerator_Test.cpp.
To select a certain distribution, hard code the identifiers and parameters in the source code of the script where it says
//
// START OF USER-DEFINED OUTPUT
//
nstates = 3;
states[0] = 0.;
states[1] = 1.;
states[2] = 0.;
mix[0] = 0.;
mix[1] = 0.;
//
// END OF USER-DEFINED OUTPUT
//
and recompile afterwards. Running the script with the -f option (the only option that is different from GetHistogram) will fit the input distribution to the simulation output.
Now, by looking at the fit
and especially the fit residuals
one can see clear systematic deviations from the input distribution which are a clear indication that W_max == 1 is not a good choice for this distribution.
Copyright (C) 2017
U. Gayer ([email protected])
O. Papst
This code is distributed under the terms of the GNU General Public License. See the COPYING file for more information.
UG would like to acknowledge the untiring effort of user Jörn Kleemann in debugging the code and making utr better in general.
[1] S. Agostinelli et al., “GEANT4 - a simulation toolkit”, Nucl. Inst. Meth. A 506.3, 250 (2003). doi:10.1016/S0168-9002(03)01368-8.
[2] J. Allison et al., “GEANT4 developments and applications”, IEEE Transactions on Nuclear Science, 53.1, 270 (2006). doi:10.1109/TNS.2006.869826.
[3] J. Allison et al., “Recent developments in GEANT4”, Nucl. Inst. Meth. A 835, 186 (2016). doi:10.1016/j.nima.2016.06.125.
[4] H. R. Weller et al., “Research opportunities at the upgraded HIγS facility”, Prog. Part. Nucl. Phys. 62.1, 257 (2009). doi:10.1016/j.ppnp.2008.07.001.
[5] B. Löher et al., “The high-efficiency γ-ray spectroscopy setup γ³ at HIγS”, Nucl. Instr. Meth. Phys. Res. A 723, 136 (2013). doi:10.1016/j.nima.2013.04.087.