Merge pull request #74 from rest-for-physics/jgalan_examples_upgrade

Adding Fe55 decay validation and new Detector Response example
rest-for-physics · Oct 7, 2022 · 1b31d60 · 1b31d60
2 parents a5f736b + 265b994
commit 1b31d60
Show file tree

Hide file tree

Showing 12 changed files with 356 additions and 2 deletions.
diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml
@@ -366,13 +366,20 @@ jobs:
         with:
           path: ${{ env.REST_PATH }}
           key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
-      - name: Run example
+      - name: Launching Gammas 
         run: |
           source ${{ env.REST_PATH }}/thisREST.sh
           cd ${{ env.REST_PATH }}/examples/restG4/11.Xrays/
           restG4 xrays.rml -o xrays_simulation.root
           restManager --c analysis.rml --f xrays_simulation.root --o xrays_simulation_analysis.root
           restRoot -b -q GetQE.C'("xrays_simulation_analysis.root")'
+      - name: Launching Fe55 source
+        run: |
+          source ${{ env.REST_PATH }}/thisREST.sh
+          cd ${{ env.REST_PATH }}/examples/restG4/11.Xrays/
+          restG4 Fe55.rml -o Fe55_simulation.root
+          restManager --c analysis.rml --f Fe55_simulation.root --o Fe55_simulation_analysis.root
+          restRoot -b -q ValidateFe55.C'("Fe55_simulation_analysis.root")'
 
   restG4-examples-12:
     name: "Example 12: Generators"
@@ -416,6 +423,28 @@ jobs:
           restG4 Neutrons.rml -o Neutrons.root -e 1
           restRoot -b -q Validate.C'("Neutrons.root")'
 
+  restG4-examples-14:
+    name: "Example 14: DetectorResponse"
+    runs-on: ubuntu-latest
+    container:
+      image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics
+    needs: [ check-installation, restG4-examples-10 ]
+    steps:
+      - uses: actions/checkout@v3
+      - name: Restore cache
+        uses: actions/cache@v3
+        id: framework-install-restG4-cache
+        with:
+          path: ${{ env.REST_PATH }}
+          key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
+      - name: Run example
+        run: |
+          source ${{ env.REST_PATH }}/thisREST.sh
+          cd ${{ env.REST_PATH }}/examples/restG4/14.DetectorResponse/
+          source launchResponse.sh
+          restManager --c analysis.rml --f RestG4_XenonNeon_30Pct_1.5bar_Drift3cm_Size6cm.root
+          restRoot -b -q ValidateResponse.C'("G4Analysis_XenonNeon_30Pct_1.5bar_Drift3cm_Size6cm.root")'
+
   # Reference version of Geant4
 
   framework-install-reference:

diff --git a/.gitignore b/.gitignore
@@ -0,0 +1 @@
+*.root
diff --git a/examples/11.Xrays/Fe55.rml b/examples/11.Xrays/Fe55.rml
@@ -0,0 +1,64 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Constructing XML-like REST metadata input file
+     File : config.rml
+     First concept author J. Galan (26-Apr-2015)
+ -->
+<!-- 
+By default REST units are mm, keV and degrees
+-->
+<restG4>
+    <globals>
+        <parameter name="mainDataPath" value=""/>
+        <parameter name="verboseLevel" value="essential"/>
+    </globals>
+    <TRestRun name="Run metadata" title="REST Metadata run info (template)">
+        <parameter name="experimentName" value="ArgonISO"/>
+        <parameter name="runType" value="simulation"/>
+        <parameter name="runNumber" value="1"/>
+        <parameter name="runTag" value="Fe55"/>
+        <parameter name="outputFileName" value="RestG4_[fRunTag]_[fRunNumber]_[fExperimentName].root"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="overwrite" value="off"/>
+        <parameter name="readOnly" value="false"/>
+    </TRestRun>
+    <TRestGeant4Metadata name="xrays" title="xrays">
+        <parameter name="gdmlFile" value="geometry/setup.gdml"/>
+        <parameter name="subEventTimeDelay" value="100" units="us"/>
+        <parameter name="seed" value="17021981"/>
+        <parameter name="nEvents" value="100000"/>
+        <parameter name="registerEmptyTracks" value="false"/>
+        <parameter name="saveAllEvents" value="true"/>
+        <generator type="point" position="(0,0,-100)" units="mm">
+            <source particle="Fe55" fullChain="on">
+                <angular type="isotropic"/>
+                <energy type="mono" energy="0keV"/>
+            </source>
+        </generator>
+        <detector>
+            <parameter name="energyRange" value="(0,30)" units="keV"/>
+            <volume name="gas" sensitive="true" maxStepSize="0.005mm"/>
+        </detector>
+    </TRestGeant4Metadata>
+    <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"/>
+
+        // Use only one EM physics list
+        <!-- EM Physics lists -->
+        <!--<physicsList name="G4EmLivermorePhysics"> </physicsList>-->
+        <!-- <physicsList name="G4EmPenelopePhysics"> </physicsList> -->
+        <physicsList name="G4EmStandardPhysics_option4"/>
+
+        <!-- Decay physics lists -->
+        <physicsList name="G4DecayPhysics"/>
+        <physicsList name="G4RadioactiveDecayPhysics"/>
+        <physicsList name="G4RadioactiveDecay"><option name="ICM" value="true"/><option name="ARM" value="true"/></physicsList>
+
+        <!-- Hadron physics lists -->
+        <physicsList name="G4HadronElasticPhysicsHP"/>
+        <physicsList name="G4IonBinaryCascadePhysics"/>
+        <physicsList name="G4HadronPhysicsQGSP_BIC_HP"/>
+        <physicsList name="G4NeutronTrackingCut"/>
+        <physicsList name="G4EmExtraPhysics"/>
+
+    </TRestGeant4PhysicsLists>
+</restG4>
diff --git a/examples/11.Xrays/README.md b/examples/11.Xrays/README.md
@@ -9,15 +9,31 @@ This example uses an x-ray gun with photons in the (0,15)keV range hitting a gas
 This example implements a simple geometry based on a cylinder full of a Xenon gas mixture.
 
 #### Testing the example
+
+There are two independent examples that exploit this detector geometry
+
+##### Example 1. Launching gammas
+
 Execute the following to launch the x-rays
 
 ```
 restG4 xrays.rml
-restManager --c analsis.rml --f RestG4_xrays_00001_ArgonISO.root
+restManager --c analysis.rml --f RestG4_xrays_00001_ArgonISO.root
 ```
 
 The script GetQE.C will produce few histograms with the calculated efficiency:
 
 ```
 restRoot -b -q GetQE.C'("Run_g4Analysis_00001_xrays.root")'
 ```
+
+##### Example 2. Fe55 isotropic source generator
+
+Execute the following to launch the Fe55 source events
+
+```
+restG4 Fe55.rml -o RestG4_Fe55_ArgonISO.root
+restManager --c analysis.rml --f RestG4_Fe55_ArgonISO.root
+```
+
+The `ValidationFe55.C` is used on the pipeline to verify that the number of events inside the peak remains correct.
diff --git a/examples/11.Xrays/ValidateFe55.C b/examples/11.Xrays/ValidateFe55.C
@@ -0,0 +1,19 @@
+Int_t ValidateFe55(std::string fName) {
+    TRestRun run(fName);
+
+    TRestAnalysisTree* aTree = run.GetAnalysisTree();
+
+    aTree->Draw("g4Ana_totalEdep");
+
+    TH1D* h = (TH1D*)aTree->GetHistogram();
+
+    Int_t peakCounts = h->Integral(h->FindFixBin(5.7), h->FindFixBin(5.9));
+
+    if (peakCounts < 1600 || peakCounts > 1750) {
+        std::cout << "Problem on Fe55 gamma peak identification." << std::endl;
+        std::cout << "Peak counts found : " << peakCounts << std::endl;
+        std::cout << "The result should be in the range ( 1600, 1750 )" << std::endl;
+    }
+
+    return 0;
+}
diff --git a/examples/14.DetectorResponse/DetectorResponse.rml b/examples/14.DetectorResponse/DetectorResponse.rml
@@ -0,0 +1,58 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Globals file contains the paths to the geometry and data directories.
+	If externally unchanged, ${REST_DATAPAH} points to current directory by default. -->
+<restG4>
+    <globals><parameter name="mainDataPath" value="."/><parameter name="verboseLevel" value="warning"/> %options are : essential silent, warning, info, debug
+	</globals>
+    <TRestRun name="QE_Xenon" title="QE detector response simulation.">
+        <parameter name="experiment" value="GasDetector"/>
+        <parameter name="runType" value="DetectorResponse"/>
+        <parameter name="runNumber" value="auto"/>
+        <parameter name="runTag" value="${GDML_GAS}_${GDML_QUENCHER_PCT}Pct_${GDML_PRESSURE}bar_Drift${GDML_DRIFT}cm_Size${GDML_DETECTOR_SIZE}cm"/>
+        <parameter name="runDescription" value="We launch photons 0 to 15 keV."/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="overwrite" value="off"/>
+        <parameter name="outputFileName" value="RestG4_[fRunTag].root"/>
+    </TRestRun>
+    <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">
+        <parameter name="gdml_file" value="geometry.gdml"/>
+        <parameter name="subEventTimeDelay" value="100" units="us"/>
+        <!-- The number of events to be simulated is now defined in TRestGeant4Metadata -->
+        <parameter name="Nevents" value="150000"/>
+        <parameter name="seed" value="137"/>
+        <parameter name="saveAllEvents" value="false"/>
+        <parameter name="printProgress" value="true"/>
+        <generator type="point" position="(0,0,-100)" units="mm">
+            <!-- postion="(0,0,100)" -->
+            <source particle="gamma">
+                <angularDist type="flux" direction="(0,0,1)"/>
+                <energyDist type="flat" range="(0,15)" units="keV"/>
+            </source>
+        </generator>
+        <detector>
+            <parameter name="energyRange" value="(0,15)" units="keV"/>
+            <volume name="target" sensitive="true" chance="1" maxStepSize="0.2mm"/>
+        </detector>
+    </TRestGeant4Metadata>
+    <TRestGeant4PhysicsLists name="default" title="First physics list implementation." verboseLevel="warning">
+        <parameter name="cutForGamma" value="1" units="um"/>
+        <parameter name="cutForElectron" value="1" units="um"/>
+        <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"/>
+        <physicsList name="G4EmLivermorePhysics"> </physicsList>
+        <physicsList name="G4DecayPhysics"> </physicsList>
+        <physicsList name="G4RadioactiveDecayPhysics"> </physicsList>
+        <physicsList name="G4RadioactiveDecay">
+            <option name="ICM" value="true"/>
+            <option name="ARM" value="true"/>
+        </physicsList>
+        <physicsList name="G4HadronElasticPhysicsHP"> </physicsList>
+        <physicsList name="G4IonBinaryCascadePhysics"> </physicsList>
+        <physicsList name="G4HadronPhysicsQGSP_BIC_HP"> </physicsList>
+        <physicsList name="G4NeutronTrackingCut"> </physicsList>
+        <physicsList name="G4EmExtraPhysics"> </physicsList>
+    </TRestGeant4PhysicsLists>
+</restG4>
diff --git a/examples/14.DetectorResponse/README.md b/examples/14.DetectorResponse/README.md
@@ -0,0 +1,53 @@
+**Date:** 30/September/2022
+
+**Author:** Javier Galan (javier.galan@unizar.es)
+
+#### Detector geometry
+
+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`.
+
+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`.
+
+#### Event generator
+
+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.
+
+#### Testing the example
+
+Execute the following to generate a dataset.
+
+```
+export GDML_PRESSURE=1.5
+export GDML_NEON_PCT=30
+export GDML_DRIFT=3
+export GDML_DETECTOR_SIZE=6
+restG4 DetectorResponse.rml 
+```
+
+or simply use the bash script
+
+```
+source launchResponse.sh
+```
+
+Then, populate the analysis tree by executing
+
+```
+restManager --c g4Analysis.rml --f Input.root --o Output.root
+```
+
+where `Input.root` is the file previously generated by `restG4`.
+
+Then open and explore the file using:
+
+```
+restRoot Output.root
+[0] new TBrowser
+```
+
+#### Results
+
+The execution of the example should produce the following result for pure xenon.
+
+![PureXenon](images/PureXenonResponse.png)
+
diff --git a/examples/14.DetectorResponse/ValidateResponse.C b/examples/14.DetectorResponse/ValidateResponse.C
@@ -0,0 +1,20 @@
+Int_t ValidateResponse(std::string fName) {
+    TRestRun run(fName);
+    TRestAnalysisTree* aTree = run.GetAnalysisTree();
+    aTree->Draw("g4Ana_totalEdep");
+    TH1D* h = (TH1D*)aTree->GetHistogram();
+
+    TRestGeant4Metadata* md = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
+
+    Double_t efficiency = h->Integral() / md->GetNumberOfEvents();
+
+    Double_t efficiencyReference = 0.91;
+    std::cout << "Overall efficiency : " << efficiency << std::endl;
+
+    if (TMath::Abs(efficiency - efficiencyReference) > 0.03) {
+        cout << "The efficiency does not match the reference value of " << efficiencyReference << endl;
+        return 1;
+    }
+
+    return 0;
+}
diff --git a/examples/14.DetectorResponse/analysis.rml b/examples/14.DetectorResponse/analysis.rml
@@ -0,0 +1,28 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!--This file is an example of REST simulation functionality. We process the output root file 
+from restG4, converting its TRestGeant4Event to TRestDetectorSignalEvent. Observables and processes's 
+internal values are saved.
+-->
+<TRestManager name="RESTManagerSim" title="Template manager to process a simulation generated by restG4.">
+    <globals><parameter name="mainDataPath" value="."/><parameter name="verboseLevel" value="warning"/> %options are : essential silent, warning, info, debug
+    </globals>
+    <TRestRun name="Process" title="Geant4 basic analysis">
+        <parameter name="experimentName" value="preserve"/>
+        <parameter name="readOnly" value="false"/>
+        <parameter name="runNumber" value="preserve"/>
+        <parameter name="runTag" value="preserve"/>
+        <parameter name="runType" value="G4Analysis"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="verboseLevel" value="1"/>
+        <parameter name="overwrite" value="off"/>
+        <parameter name="outputFileName" value="[fRunType]_[fRunTag].root"/>
+    </TRestRun>
+    <TRestProcessRunner name="TemplateEventProcess" verboseLevel="info"><parameter name="eventsToProcess" value="0"/><parameter name="inputEventStorage" value="OFF"/><parameter name="outputEventStorage" value="OFF"/>
+
+        // observables = all will add all NON `custom` observables
+		<addProcess type="TRestGeant4AnalysisProcess" name="g4Ana" observable="all"/>
+
+    </TRestProcessRunner>
+    <addTask type="processEvents" value="ON"/>
+</TRestManager>
diff --git a/examples/14.DetectorResponse/geometry.gdml b/examples/14.DetectorResponse/geometry.gdml
@@ -0,0 +1,56 @@
+<?xml version="1.0" encoding="utf-8" standalone="no" ?>
+<!-- ##VERSION REST ${GDML_GAS} (${GDML_QUENCHER_PCT}Pct) - ${GDML_PRESSURE}bar - Drift:${GDML_DRIFT}cm - Detector size:${GDML_DETECTOR_SIZE}cm ## -->
+
+<!DOCTYPE gdml [
+        <!ENTITY geometry SYSTEM "geometry.gdml">
+		<!ENTITY materials SYSTEM "https://rest-for-physics.github.io/materials/rest.xml">
+        ]>
+
+<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+      xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
+
+    <define>
+		<!-- All lenghts should be in mm -->
+        <constant name="world_size" value="20000"/>
+
+		<variable name="targetGasDensity" value="(100-${GDML_QUENCHER_PCT})/100. * ${GDML_PRESSURE} * ${GDML_TARGET_DENSITY} * (273.15+15) / 273.15 / 1.01325"/> 
+
+		<variable name="quencherDensity" value="${GDML_QUENCHER_PCT} / 100. * ${GDML_PRESSURE} * ${GDML_QUENCHER_DENSITY}  * (273.15+15) / 273.15 / 1.01325"/>
+		<variable name="quencherFraction" value="${GDML_QUENCHER_PCT} / 100."/>
+
+		<!-- What matters is just the density (this is just additional info) -->
+        <variable name="gasTemperature" value="273.15+15"/>
+		<variable name="gasPressure" value="${GDML_PRESSURE}"/>
+    </define>
+
+    &materials;
+
+	<solids>
+		<box name="WorldSolid" x="world_size" y="world_size" z="world_size" lunit="mm"/>
+		<box name="targetSolid" startphi="0" deltaphi="360" x="${GDML_DETECTOR_SIZE}" y="${GDML_DETECTOR_SIZE}" z="${GDML_DRIFT}" lunit="cm"/>
+	</solids>
+
+	<structure>
+
+	<volume name="targetVolume">
+		<materialref ref="${GDML_GAS}"/>
+		<solidref ref="targetSolid"/>
+	</volume>
+
+	<volume name="World">
+		<materialref ref="Vacuum"/>
+		<solidref ref="WorldSolid"/>
+
+		<physvol name="target">
+			<volumeref ref="targetVolume"/>
+			<position name="targetPos" unit="mm" x="0" y="0" z="0"/>
+		</physvol>
+	</volume>
+
+	</structure>
+
+    <setup name="Default" version="1.0">
+        <world ref="World"/>
+    </setup>
+
+</gdml>
diff --git a/examples/14.DetectorResponse/images/PureXenonResponse.png b/examples/14.DetectorResponse/images/PureXenonResponse.png