Skip to content

Commit

Permalink
Merge pull request #28 from rest-for-physics/alphas_example
Browse files Browse the repository at this point in the history
Alphas example
  • Loading branch information
jgalan authored Jan 29, 2022
2 parents e45a940 + f76b5e2 commit c86aa2f
Show file tree
Hide file tree
Showing 14 changed files with 453 additions and 9 deletions.
35 changes: 28 additions & 7 deletions .gitlab-ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ loadRESTLibs:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/01.NLDBD/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/01.NLDBD/
- restG4 NLDBD.rml
- geant4-config --version
- restRoot -b -q Validate.C'("Run00001_NLDBD_Test.root")'
Expand All @@ -92,7 +92,7 @@ loadRESTLibs:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/03.Fluorescence/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/03.Fluorescence/
- restG4 gamma.rml
- restManager --c g4Analysis.rml --f Run00111_Gamma_Fluorescence.root
- geant4-config --version
Expand All @@ -102,7 +102,7 @@ loadRESTLibs:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/04.MuonScan/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/04.MuonScan/
- restG4 Muon.rml
- geant4-config --version
- restRoot -b -q Validate.C'("Run00111_restG4_MuonsFromPoint_rest.root")'
Expand All @@ -113,23 +113,23 @@ loadRESTLibs:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/05.PandaXIII/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/05.PandaXIII/
- restG4 Xe136bb0n.rml
- restRoot -b -q ValidateG4.C'("Xe136bb0n_n2E06.root")'

06_IonRecoils:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/06.IonRecoils/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/06.IonRecoils/
- restG4 recoils.rml
- restRoot -b -q Validate.C'("Run00001_F20_Recoils.root")'

07_FullChainDecay:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/07.FullChainDecay/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/07.FullChainDecay/
- restG4 fullChain.rml
- restG4 singleDecay.rml
- restRoot -b -q Validate.C'("Run00001_U238_FullChainDecay.root", 14)'
Expand All @@ -138,11 +138,32 @@ loadRESTLibs:
- restG4 singleDecay.rml
- restRoot -b -q Validate.C'("Run00002_Be7_SingleChainDecay.root", 1)'

08_Alphas:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/restG4/08.Alphas/
- mkdir data
- export REST_ENERGY=5
- export REST_FOIL=1
- restG4 alphas.rml
- export REST_ENERGY=5
- export REST_FOIL=5
- restG4 alphas.rml
- export REST_ENERGY=1
- export REST_FOIL=5
- restG4 alphas.rml
- restManager --c analysis.rml --f data/Run_1MeV_5um.root
- restManager --c analysis.rml --f data/Run_5MeV_1um.root
- restManager --c analysis.rml --f data/Run_5MeV_5um.root
- restManager --c plots.rml --f "data/*g4Ana*root" --batch
- restRoot -b -q Validate.C


09_Pb210_Shielding:
type: examples
script:
- . ${CI_PROJECT_DIR}/install/thisREST.sh
- cd ${CI_PROJECT_DIR}/install/examples/09.Pb210_Shield/
- cd ${CI_PROJECT_DIR}/install/examples/restG4/09.Pb210_Shield/
- restG4 Pb210.rml
- restRoot -b -q Validate.C'("Run00001_Pb210_Shielding.root")'
4 changes: 2 additions & 2 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,8 @@ target_link_libraries(${PROJECT_NAME} ${LINK_LIBRARIES})
target_include_directories(${PROJECT_NAME} SYSTEM PUBLIC ${INCLUDE_DIRS})

# Copy scripts to the build directory
install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/examples
DESTINATION .
install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/examples/
DESTINATION ./examples/restG4/
COMPONENT install
)

Expand Down
81 changes: 81 additions & 0 deletions examples/08.Alphas/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
**Date:** 25/January/2022

**Author:** Javier Galan (javier.galan@unizar.es)

This example tests the generation of alphas with a well defined energy from a thin foil immersed in a gas volume.

#### Detector geometry

The geometry consists of a 10x10x10cm3 box made of Argon gas. In the middle of the box a small copper plate is placed. The thickness of the plate in um can be controlled with an enviromental system variable, REST_FOIL.

#### Event generator
The event generator will produce mono-energetic alphas produced in the copper sample volume (physical volume named source). The energy of the alphas can be controled using the REST_ENERGY system environment variable expressed in MeV.

#### Testing the example

Create a "data" directory if not existing already

```
mkdir data
```

Execute the following to launch the generation of 10,000 alphas.

```
restG4 alphas.rml
```

Execute the following to launch the processing of the generated events and produce few observables inside the analysis tree.

```
restG4 --c analysis.rml --f data/FileGenerated.root
```

See also `../../.gitlab-ci.yml` and `Validation.C` for validation routines executed in the pipeline. The following plot showing the energy spectrum is generated using the pipeline.


![plot](plots.png)

Commands to reproduce the figure shown follow:

```
mkdir data
export REST_ENERGY=5
export REST_FOIL=1
restG4 alphas.rml
export REST_FOIL=5
restG4 alphas.rml
export REST_ENERGY=1
restG4 alphas.rml
restManager --c analysis.rml --f Run_1MeV_5um.root
restManager --c analysis.rml --f Run_5MeV_1um.root
restManager --c analysis.rml --f Run_5MeV_5um.root
restManager --c plots.rml --f "data/*g4Ana*root"
```


A minimal ROOT shell commands to visualize the detector geometry is:

```
export REST_FOIL=5
restRoot
root [0] TRestGDMLParser *g = new TRestGDMLParser()
root [1] TGeoManager *geo = g->GetGeoManager("geometry/setup.gdml")
root [2] geo->GetTopVolume()->Draw("ogl")
```

![geometry](geometry.png)

Minimal ROOT shell commands to visualize events

```
restRootMacros
root [0] REST_Geant4_ViewEvent("data/Run_5MeV_1um.root")
root [1] TGeoManager *geo = g->GetGeoManager("geometry/setup.gdml")
root [2] geo->GetTopVolume()->Draw("ogl")
```

![event](event.png)
61 changes: 61 additions & 0 deletions examples/08.Alphas/Validate.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@

Int_t Validate() {
TFile* f = new TFile("plots.root");

TH1F* h1 = (TH1F*)f->Get("1MeV_5um");
if (h1->Integral() != 1578) {
cout << "Wrong number of alphas produced at 1MeV_5um!" << endl;
cout << "Histogram contains : " << h1->Integral() << endl;
return 10;
}

TH1F* h2 = (TH1F*)f->Get("5MeV_5um");
if (h2->Integral() != 7182) {
cout << "Wrong number of alphas produced at 5MeV_5um!" << endl;
cout << "Histogram contains : " << h2->Integral() << endl;
return 20;
}

TH1F* h3 = (TH1F*)f->Get("5MeV_1um");
if (h3->Integral() != 9421) {
cout << "Wrong number of alphas produced at 5MeV_1um!" << endl;
cout << "Histogram contains : " << h3->Integral() << endl;
return 30;
}

TRestRun* run = new TRestRun("data/Run_g4Analysis_1MeV_5um.root");
TRestAnalysisTree* aT = run->GetAnalysisTree();

run->GetEntry(100);

Int_t thetaInt = (int)(1000 * aT->GetObservableValue<Double_t>("g4Ana_thetaPrimary"));
if (thetaInt != 515) {
cout << "Wrong theta angle value for entry 100!" << endl;
cout << "Theta value is : " << thetaInt << endl;
return 60;
}

Int_t phiInt = (int)(1000 * aT->GetObservableValue<Double_t>("g4Ana_phiPrimary"));
if (phiInt != 1539) {
cout << "Wrong phi angle value for entry 100!" << endl;
cout << "Phi value is : " << phiInt << endl;
return 60;
}

Double_t theta = aT->GetObservableAverage("g4Ana_thetaPrimary");
if ((int)(1000 * theta) != 1589) {
cout << "Wrong theta angle average!" << endl;
cout << "Theta angle average : " << theta << " while it should be : " << 1.589 << endl;
return 80;
}

Double_t phi = aT->GetObservableAverage("g4Ana_phiPrimary");
if ((int)(1000 * phi) != 27) {
cout << "Wrong phi angle average!" << endl;
cout << "Phi angle average : " << phi << " while it should be : " << 0.027 << endl;
return 90;
}

cout << "All tests passed! [\033[32mOK\033[0m]\n";
return 0;
}
90 changes: 90 additions & 0 deletions examples/08.Alphas/alphas.rml
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
<?xml version="4.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="Alphas"/>
<parameter name="runType" value="simulation"/>
<parameter name="runNumber" value="1"/>
<parameter name="runTag" value="${REST_ENERGY}MeV_${REST_FOIL}um"/>
<parameter name="outputFileName" value="data/Run_[fRunTag].root"/>
<parameter name="runDescription" value=""/>
<parameter name="user" value="${USER}"/>
<parameter name="overwrite" value="off"/>
<parameter name="readOnly" value="false"/>
</TRestRun>

<TRestGeant4Metadata name="alphas" title="Alphas">

<parameter name="gdml_file" value="geometry/setup.gdml"/>
<parameter name="subEventTimeDelay" value="100" units="us"/>

<parameter name="seed" value="17021981"/>

<!-- The number of events to be simulated is now defined in TRestGeant4Metadata -->
<parameter name="Nevents" value="10000"/>
<parameter name="registerEmptyTracks" value="false"/>
<parameter name="saveAllEvents" value="false"/>

<generator type="volume" from="source" >
<source particle="alpha">
<angularDist type="isotropic"/>
<energyDist type="mono" energy="${REST_ENERGY}" units="MeV" />
</source>
</generator>

<storage sensitiveVolume="gas">
<parameter name="energyRange" value="(0,2*${REST_ENERGY})" units="MeV"/>
<activeVolume name="gas" chance="1" maxStepSize="0.5mm"/>
</storage>
</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"></physicsList>

<!-- Decay physics lists -->
<physicsList name="G4DecayPhysics"></physicsList>
<physicsList name="G4RadioactiveDecayPhysics"></physicsList>
<physicsList name="G4RadioactiveDecay">
<option name="ICM" value="true"/>
<option name="ARM" value="true"/>
</physicsList>

<!-- Hadron physics lists -->
<physicsList name="G4HadronElasticPhysicsHP"></physicsList>
<physicsList name="G4IonBinaryCascadePhysics"></physicsList>
<physicsList name="G4HadronPhysicsQGSP_BIC_HP"></physicsList>
<physicsList name="G4NeutronTrackingCut"></physicsList>
<physicsList name="G4EmExtraPhysics"></physicsList>

</TRestGeant4PhysicsLists>

</restG4>
53 changes: 53 additions & 0 deletions examples/08.Alphas/analysis.rml
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
<?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>
<variable name="VERSION" value="1.0" overwrite="false"/>
<variable name="EXPERIMENT" value="TREXDM" overwrite="false"/>
<parameter name="mainDataPath" value="."/>
<parameter name="verboseLevel" value="warning"/>
</globals>

<TRestRun name="Process" title="${EXPERIMENT} Simulations. Version ${VERSION}.">
<parameter name="experimentName" value="${EXPERIMENT}"/>
<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="data/Run_[fRunType]_[fRunTag].root"/>
</TRestRun>

<TRestProcessRunner name="TemplateEventProcess" verboseLevel="info">
<parameter name="eventsToProcess" value="0"/>
<parameter name="threadNumber" value="2"/>

<parameter name="inputAnalysisStorage" value="on"/>
<parameter name="inputEventStorage" value="off"/>
<parameter name="outputEventStorage" value="off"/>

// observable = all will add all NON `custom` observables
<addProcess type="TRestGeant4AnalysisProcess" name="g4Ana" value="ON">
<observable name="thetaPrimary" />
<observable name="phiPrimary" />
<observable name="xOriginPrimary" />
<observable name="yOriginPrimary" />
<observable name="zOriginPrimary" />
<observable name="totalEdep" />
</addProcess>

</TRestProcessRunner>

<addTask type="processEvents" value="ON"/>

</TRestManager>

Binary file added examples/08.Alphas/event.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/08.Alphas/geometry.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading

0 comments on commit c86aa2f

Please sign in to comment.