Skip to content

Commit 1b31d60

Browse files
authored
Merge pull request #74 from rest-for-physics/jgalan_examples_upgrade
Adding Fe55 decay validation and new Detector Response example
2 parents a5f736b + 265b994 commit 1b31d60

File tree

12 files changed

+356
-2
lines changed

12 files changed

+356
-2
lines changed

.github/workflows/validation.yml

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -366,13 +366,20 @@ jobs:
366366
with:
367367
path: ${{ env.REST_PATH }}
368368
key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
369-
- name: Run example
369+
- name: Launching Gammas
370370
run: |
371371
source ${{ env.REST_PATH }}/thisREST.sh
372372
cd ${{ env.REST_PATH }}/examples/restG4/11.Xrays/
373373
restG4 xrays.rml -o xrays_simulation.root
374374
restManager --c analysis.rml --f xrays_simulation.root --o xrays_simulation_analysis.root
375375
restRoot -b -q GetQE.C'("xrays_simulation_analysis.root")'
376+
- name: Launching Fe55 source
377+
run: |
378+
source ${{ env.REST_PATH }}/thisREST.sh
379+
cd ${{ env.REST_PATH }}/examples/restG4/11.Xrays/
380+
restG4 Fe55.rml -o Fe55_simulation.root
381+
restManager --c analysis.rml --f Fe55_simulation.root --o Fe55_simulation_analysis.root
382+
restRoot -b -q ValidateFe55.C'("Fe55_simulation_analysis.root")'
376383
377384
restG4-examples-12:
378385
name: "Example 12: Generators"
@@ -416,6 +423,28 @@ jobs:
416423
restG4 Neutrons.rml -o Neutrons.root -e 1
417424
restRoot -b -q Validate.C'("Neutrons.root")'
418425
426+
restG4-examples-14:
427+
name: "Example 14: DetectorResponse"
428+
runs-on: ubuntu-latest
429+
container:
430+
image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics
431+
needs: [ check-installation, restG4-examples-10 ]
432+
steps:
433+
- uses: actions/checkout@v3
434+
- name: Restore cache
435+
uses: actions/cache@v3
436+
id: framework-install-restG4-cache
437+
with:
438+
path: ${{ env.REST_PATH }}
439+
key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
440+
- name: Run example
441+
run: |
442+
source ${{ env.REST_PATH }}/thisREST.sh
443+
cd ${{ env.REST_PATH }}/examples/restG4/14.DetectorResponse/
444+
source launchResponse.sh
445+
restManager --c analysis.rml --f RestG4_XenonNeon_30Pct_1.5bar_Drift3cm_Size6cm.root
446+
restRoot -b -q ValidateResponse.C'("G4Analysis_XenonNeon_30Pct_1.5bar_Drift3cm_Size6cm.root")'
447+
419448
# Reference version of Geant4
420449

421450
framework-install-reference:

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
*.root

examples/11.Xrays/Fe55.rml

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
2+
<!-- Constructing XML-like REST metadata input file
3+
File : config.rml
4+
First concept author J. Galan (26-Apr-2015)
5+
-->
6+
<!--
7+
By default REST units are mm, keV and degrees
8+
-->
9+
<restG4>
10+
<globals>
11+
<parameter name="mainDataPath" value=""/>
12+
<parameter name="verboseLevel" value="essential"/>
13+
</globals>
14+
<TRestRun name="Run metadata" title="REST Metadata run info (template)">
15+
<parameter name="experimentName" value="ArgonISO"/>
16+
<parameter name="runType" value="simulation"/>
17+
<parameter name="runNumber" value="1"/>
18+
<parameter name="runTag" value="Fe55"/>
19+
<parameter name="outputFileName" value="RestG4_[fRunTag]_[fRunNumber]_[fExperimentName].root"/>
20+
<parameter name="runDescription" value=""/>
21+
<parameter name="user" value="${USER}"/>
22+
<parameter name="overwrite" value="off"/>
23+
<parameter name="readOnly" value="false"/>
24+
</TRestRun>
25+
<TRestGeant4Metadata name="xrays" title="xrays">
26+
<parameter name="gdmlFile" value="geometry/setup.gdml"/>
27+
<parameter name="subEventTimeDelay" value="100" units="us"/>
28+
<parameter name="seed" value="17021981"/>
29+
<parameter name="nEvents" value="100000"/>
30+
<parameter name="registerEmptyTracks" value="false"/>
31+
<parameter name="saveAllEvents" value="true"/>
32+
<generator type="point" position="(0,0,-100)" units="mm">
33+
<source particle="Fe55" fullChain="on">
34+
<angular type="isotropic"/>
35+
<energy type="mono" energy="0keV"/>
36+
</source>
37+
</generator>
38+
<detector>
39+
<parameter name="energyRange" value="(0,30)" units="keV"/>
40+
<volume name="gas" sensitive="true" maxStepSize="0.005mm"/>
41+
</detector>
42+
</TRestGeant4Metadata>
43+
<TRestGeant4PhysicsLists name="default" title="First physics list implementation."><parameter name="cutForGamma" value="0.01" units="mm"/><parameter name="cutForElectron" value="2" units="mm"/><parameter name="cutForPositron" value="1" units="mm"/><parameter name="cutForMuon" value="1" units="mm"/><parameter name="cutForNeutron" value="1" units="mm"/><parameter name="minEnergyRangeProductionCuts" value="1" units="keV"/><parameter name="maxEnergyRangeProductionCuts" value="1" units="GeV"/>
44+
45+
// Use only one EM physics list
46+
<!-- EM Physics lists -->
47+
<!--<physicsList name="G4EmLivermorePhysics"> </physicsList>-->
48+
<!-- <physicsList name="G4EmPenelopePhysics"> </physicsList> -->
49+
<physicsList name="G4EmStandardPhysics_option4"/>
50+
51+
<!-- Decay physics lists -->
52+
<physicsList name="G4DecayPhysics"/>
53+
<physicsList name="G4RadioactiveDecayPhysics"/>
54+
<physicsList name="G4RadioactiveDecay"><option name="ICM" value="true"/><option name="ARM" value="true"/></physicsList>
55+
56+
<!-- Hadron physics lists -->
57+
<physicsList name="G4HadronElasticPhysicsHP"/>
58+
<physicsList name="G4IonBinaryCascadePhysics"/>
59+
<physicsList name="G4HadronPhysicsQGSP_BIC_HP"/>
60+
<physicsList name="G4NeutronTrackingCut"/>
61+
<physicsList name="G4EmExtraPhysics"/>
62+
63+
</TRestGeant4PhysicsLists>
64+
</restG4>

examples/11.Xrays/README.md

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,15 +9,31 @@ This example uses an x-ray gun with photons in the (0,15)keV range hitting a gas
99
This example implements a simple geometry based on a cylinder full of a Xenon gas mixture.
1010

1111
#### Testing the example
12+
13+
There are two independent examples that exploit this detector geometry
14+
15+
##### Example 1. Launching gammas
16+
1217
Execute the following to launch the x-rays
1318

1419
```
1520
restG4 xrays.rml
16-
restManager --c analsis.rml --f RestG4_xrays_00001_ArgonISO.root
21+
restManager --c analysis.rml --f RestG4_xrays_00001_ArgonISO.root
1722
```
1823

1924
The script GetQE.C will produce few histograms with the calculated efficiency:
2025

2126
```
2227
restRoot -b -q GetQE.C'("Run_g4Analysis_00001_xrays.root")'
2328
```
29+
30+
##### Example 2. Fe55 isotropic source generator
31+
32+
Execute the following to launch the Fe55 source events
33+
34+
```
35+
restG4 Fe55.rml -o RestG4_Fe55_ArgonISO.root
36+
restManager --c analysis.rml --f RestG4_Fe55_ArgonISO.root
37+
```
38+
39+
The `ValidationFe55.C` is used on the pipeline to verify that the number of events inside the peak remains correct.

examples/11.Xrays/ValidateFe55.C

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
Int_t ValidateFe55(std::string fName) {
2+
TRestRun run(fName);
3+
4+
TRestAnalysisTree* aTree = run.GetAnalysisTree();
5+
6+
aTree->Draw("g4Ana_totalEdep");
7+
8+
TH1D* h = (TH1D*)aTree->GetHistogram();
9+
10+
Int_t peakCounts = h->Integral(h->FindFixBin(5.7), h->FindFixBin(5.9));
11+
12+
if (peakCounts < 1600 || peakCounts > 1750) {
13+
std::cout << "Problem on Fe55 gamma peak identification." << std::endl;
14+
std::cout << "Peak counts found : " << peakCounts << std::endl;
15+
std::cout << "The result should be in the range ( 1600, 1750 )" << std::endl;
16+
}
17+
18+
return 0;
19+
}
Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
2+
<!-- Globals file contains the paths to the geometry and data directories.
3+
If externally unchanged, ${REST_DATAPAH} points to current directory by default. -->
4+
<restG4>
5+
<globals><parameter name="mainDataPath" value="."/><parameter name="verboseLevel" value="warning"/> %options are : essential silent, warning, info, debug
6+
</globals>
7+
<TRestRun name="QE_Xenon" title="QE detector response simulation.">
8+
<parameter name="experiment" value="GasDetector"/>
9+
<parameter name="runType" value="DetectorResponse"/>
10+
<parameter name="runNumber" value="auto"/>
11+
<parameter name="runTag" value="${GDML_GAS}_${GDML_QUENCHER_PCT}Pct_${GDML_PRESSURE}bar_Drift${GDML_DRIFT}cm_Size${GDML_DETECTOR_SIZE}cm"/>
12+
<parameter name="runDescription" value="We launch photons 0 to 15 keV."/>
13+
<parameter name="user" value="${USER}"/>
14+
<parameter name="overwrite" value="off"/>
15+
<parameter name="outputFileName" value="RestG4_[fRunTag].root"/>
16+
</TRestRun>
17+
<TRestGeant4Metadata name="${GDML_GAS}_${GDML_PRESSURE}bar_Drift${GDML_DRIFT}cm_Size${GDML_DETECTOR_SIZE}cm" title="${GDML_GAS} at ${GDML_PRESSURE}bar. Detector size: ${GDML_DETECTOR_SIZE}x${GDML_DETECTOR_SIZE}cm^2 and ${GDML_DRIFT}cm drift" verboseLevel="warning">
18+
<parameter name="gdml_file" value="geometry.gdml"/>
19+
<parameter name="subEventTimeDelay" value="100" units="us"/>
20+
<!-- The number of events to be simulated is now defined in TRestGeant4Metadata -->
21+
<parameter name="Nevents" value="150000"/>
22+
<parameter name="seed" value="137"/>
23+
<parameter name="saveAllEvents" value="false"/>
24+
<parameter name="printProgress" value="true"/>
25+
<generator type="point" position="(0,0,-100)" units="mm">
26+
<!-- postion="(0,0,100)" -->
27+
<source particle="gamma">
28+
<angularDist type="flux" direction="(0,0,1)"/>
29+
<energyDist type="flat" range="(0,15)" units="keV"/>
30+
</source>
31+
</generator>
32+
<detector>
33+
<parameter name="energyRange" value="(0,15)" units="keV"/>
34+
<volume name="target" sensitive="true" chance="1" maxStepSize="0.2mm"/>
35+
</detector>
36+
</TRestGeant4Metadata>
37+
<TRestGeant4PhysicsLists name="default" title="First physics list implementation." verboseLevel="warning">
38+
<parameter name="cutForGamma" value="1" units="um"/>
39+
<parameter name="cutForElectron" value="1" units="um"/>
40+
<parameter name="cutForPositron" value="1" units="mm"/>
41+
<parameter name="cutForMuon" value="1" units="mm"/>
42+
<parameter name="cutForNeutron" value="1" units="mm"/>
43+
<parameter name="minEnergyRangeProductionCuts" value="1" units="keV"/>
44+
<parameter name="maxEnergyRangeProductionCuts" value="1" units="GeV"/>
45+
<physicsList name="G4EmLivermorePhysics"> </physicsList>
46+
<physicsList name="G4DecayPhysics"> </physicsList>
47+
<physicsList name="G4RadioactiveDecayPhysics"> </physicsList>
48+
<physicsList name="G4RadioactiveDecay">
49+
<option name="ICM" value="true"/>
50+
<option name="ARM" value="true"/>
51+
</physicsList>
52+
<physicsList name="G4HadronElasticPhysicsHP"> </physicsList>
53+
<physicsList name="G4IonBinaryCascadePhysics"> </physicsList>
54+
<physicsList name="G4HadronPhysicsQGSP_BIC_HP"> </physicsList>
55+
<physicsList name="G4NeutronTrackingCut"> </physicsList>
56+
<physicsList name="G4EmExtraPhysics"> </physicsList>
57+
</TRestGeant4PhysicsLists>
58+
</restG4>
Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
**Date:** 30/September/2022
2+
3+
**Author:** Javier Galan ([email protected])
4+
5+
#### Detector geometry
6+
7+
This example implements a simple geometry composed by a `detector` box representing a gaseous TPC active volume. The detector has a squared shaped defined by the system variable `GDML_DETECTOR_SIZE` in cm, and a drift distance (defined on the z coordinate) of size `GDML_DRIFT`.
8+
9+
The geometry defines a gas mixture made of Neon and Xenon, the pressure and composition of this mixture can be controlled using `GDML_PRESSURE` in bar, and `GDML_NEON_PCT`.
10+
11+
#### Event generator
12+
13+
The event generator will launch photons parallel to the z-axis towards the center of the detector. The energy of those photons will be homogeneously distributed between 0 and 15 keV.
14+
15+
#### Testing the example
16+
17+
Execute the following to generate a dataset.
18+
19+
```
20+
export GDML_PRESSURE=1.5
21+
export GDML_NEON_PCT=30
22+
export GDML_DRIFT=3
23+
export GDML_DETECTOR_SIZE=6
24+
restG4 DetectorResponse.rml
25+
```
26+
27+
or simply use the bash script
28+
29+
```
30+
source launchResponse.sh
31+
```
32+
33+
Then, populate the analysis tree by executing
34+
35+
```
36+
restManager --c g4Analysis.rml --f Input.root --o Output.root
37+
```
38+
39+
where `Input.root` is the file previously generated by `restG4`.
40+
41+
Then open and explore the file using:
42+
43+
```
44+
restRoot Output.root
45+
[0] new TBrowser
46+
```
47+
48+
#### Results
49+
50+
The execution of the example should produce the following result for pure xenon.
51+
52+
![PureXenon](images/PureXenonResponse.png)
53+
Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
Int_t ValidateResponse(std::string fName) {
2+
TRestRun run(fName);
3+
TRestAnalysisTree* aTree = run.GetAnalysisTree();
4+
aTree->Draw("g4Ana_totalEdep");
5+
TH1D* h = (TH1D*)aTree->GetHistogram();
6+
7+
TRestGeant4Metadata* md = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
8+
9+
Double_t efficiency = h->Integral() / md->GetNumberOfEvents();
10+
11+
Double_t efficiencyReference = 0.91;
12+
std::cout << "Overall efficiency : " << efficiency << std::endl;
13+
14+
if (TMath::Abs(efficiency - efficiencyReference) > 0.03) {
15+
cout << "The efficiency does not match the reference value of " << efficiencyReference << endl;
16+
return 1;
17+
}
18+
19+
return 0;
20+
}
Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
2+
<!--This file is an example of REST simulation functionality. We process the output root file
3+
from restG4, converting its TRestGeant4Event to TRestDetectorSignalEvent. Observables and processes's
4+
internal values are saved.
5+
-->
6+
<TRestManager name="RESTManagerSim" title="Template manager to process a simulation generated by restG4.">
7+
<globals><parameter name="mainDataPath" value="."/><parameter name="verboseLevel" value="warning"/> %options are : essential silent, warning, info, debug
8+
</globals>
9+
<TRestRun name="Process" title="Geant4 basic analysis">
10+
<parameter name="experimentName" value="preserve"/>
11+
<parameter name="readOnly" value="false"/>
12+
<parameter name="runNumber" value="preserve"/>
13+
<parameter name="runTag" value="preserve"/>
14+
<parameter name="runType" value="G4Analysis"/>
15+
<parameter name="runDescription" value=""/>
16+
<parameter name="user" value="${USER}"/>
17+
<parameter name="verboseLevel" value="1"/>
18+
<parameter name="overwrite" value="off"/>
19+
<parameter name="outputFileName" value="[fRunType]_[fRunTag].root"/>
20+
</TRestRun>
21+
<TRestProcessRunner name="TemplateEventProcess" verboseLevel="info"><parameter name="eventsToProcess" value="0"/><parameter name="inputEventStorage" value="OFF"/><parameter name="outputEventStorage" value="OFF"/>
22+
23+
// observables = all will add all NON `custom` observables
24+
<addProcess type="TRestGeant4AnalysisProcess" name="g4Ana" observable="all"/>
25+
26+
</TRestProcessRunner>
27+
<addTask type="processEvents" value="ON"/>
28+
</TRestManager>
Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
<?xml version="1.0" encoding="utf-8" standalone="no" ?>
2+
<!-- ##VERSION REST ${GDML_GAS} (${GDML_QUENCHER_PCT}Pct) - ${GDML_PRESSURE}bar - Drift:${GDML_DRIFT}cm - Detector size:${GDML_DETECTOR_SIZE}cm ## -->
3+
4+
<!DOCTYPE gdml [
5+
<!ENTITY geometry SYSTEM "geometry.gdml">
6+
<!ENTITY materials SYSTEM "https://rest-for-physics.github.io/materials/rest.xml">
7+
]>
8+
9+
<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
10+
xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
11+
12+
<define>
13+
<!-- All lenghts should be in mm -->
14+
<constant name="world_size" value="20000"/>
15+
16+
<variable name="targetGasDensity" value="(100-${GDML_QUENCHER_PCT})/100. * ${GDML_PRESSURE} * ${GDML_TARGET_DENSITY} * (273.15+15) / 273.15 / 1.01325"/>
17+
18+
<variable name="quencherDensity" value="${GDML_QUENCHER_PCT} / 100. * ${GDML_PRESSURE} * ${GDML_QUENCHER_DENSITY} * (273.15+15) / 273.15 / 1.01325"/>
19+
<variable name="quencherFraction" value="${GDML_QUENCHER_PCT} / 100."/>
20+
21+
<!-- What matters is just the density (this is just additional info) -->
22+
<variable name="gasTemperature" value="273.15+15"/>
23+
<variable name="gasPressure" value="${GDML_PRESSURE}"/>
24+
</define>
25+
26+
&materials;
27+
28+
<solids>
29+
<box name="WorldSolid" x="world_size" y="world_size" z="world_size" lunit="mm"/>
30+
<box name="targetSolid" startphi="0" deltaphi="360" x="${GDML_DETECTOR_SIZE}" y="${GDML_DETECTOR_SIZE}" z="${GDML_DRIFT}" lunit="cm"/>
31+
</solids>
32+
33+
<structure>
34+
35+
<volume name="targetVolume">
36+
<materialref ref="${GDML_GAS}"/>
37+
<solidref ref="targetSolid"/>
38+
</volume>
39+
40+
<volume name="World">
41+
<materialref ref="Vacuum"/>
42+
<solidref ref="WorldSolid"/>
43+
44+
<physvol name="target">
45+
<volumeref ref="targetVolume"/>
46+
<position name="targetPos" unit="mm" x="0" y="0" z="0"/>
47+
</physvol>
48+
</volume>
49+
50+
</structure>
51+
52+
<setup name="Default" version="1.0">
53+
<world ref="World"/>
54+
</setup>
55+
56+
</gdml>

0 commit comments

Comments
 (0)