New example to generate decay spectra as a tree

examples/15.DecaySpectra/GenerateTree.C

@@ -0,0 +1,126 @@
+#include <filesystem>
+#include <memory>
+
+using namespace std;
+
+void GenerateTree(const char* filenameOrDirectoryWithFiles) {
+    gSystem->Load("libRestFramework.so");
+    gSystem->Load("libRestGeant4.so");
+
+    vector<string> filenames;
+    if (filesystem::is_directory(filenameOrDirectoryWithFiles)) {
+        for (const auto& entry : filesystem::directory_iterator(filenameOrDirectoryWithFiles)) {
+            filenames.push_back(entry.path());
+            cout << "Adding file: " << entry.path() << endl;
+        }
+    } else {
+        filenames.push_back(filenameOrDirectoryWithFiles);
+    }
+
+    TFile* outputFile = nullptr;
+    const map<string, string> particleToString = {
+        {"gamma", "gammas"},
+        /*
+        {"e-",      "electrons"},
+        {"e+",      "positrons"},
+        {"proton",  "protons"},
+        {"neutron", "neutrons"},
+        {"alpha",   "alphas"},
+         */
+    };
+    vector<string> particlesToRecord;
+    for (const auto& particle : particleToString) {
+        particlesToRecord.push_back(particle.first);
+    }
+    map<string, TTree*> trees;
+    map<string, double> energies;
+
+    bool first = true;
+    for (const auto& filename : filenames) {
+        cout << "Processing file: " << filename << endl;
+
+        TRestRun run(filename);
+        cout << "Number of entries: " << run.GetEntries() << endl;
+
+        TRestGeant4Event* event = run.GetInputEvent<TRestGeant4Event>();
+        run.GetEntry(0);
+
+        if (first) {
+            const TString primaryName = event->GetPrimaryEventParticleName();
+            // Create a new tree to store gammas from the primary particle
+            outputFile =
+                new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE");
+
+            // create a new tree for each particle
+            for (const auto& particle : particlesToRecord) {
+                // we cannot use the particle name (e.g. gamma) as a tree name, because it is a reserved word
+                trees[particle] =
+                    new TTree(particleToString.at(particle).c_str(), particleToString.at(particle).c_str());
+            }
+
+            // create a branch for each particle
+            for (const auto& particle : particlesToRecord) {
+                energies[particle] = 0;
+                trees[particle]->Branch("energy", &energies[particle], "energy/D");
+            }
+            first = false;
+        }
+
+        const string detectorName = "detector";
+        cout << "Detector name: " << detectorName << endl;
+        cout << "Starting loop over events..." << endl;
+        for (int n = 0; n < run.GetEntries(); n++) {
+            // print the progress every 1%
+            if (n % (run.GetEntries() / 100) == 0) {
+                cout << "Progress: " << n * 100 / run.GetEntries() << "%" << endl;
+            }
+            run.GetEntry(n);
+            // iterate over all particles in the event
+            for (const auto& track : event->GetTracks()) {
+                // if particle not in the list of particles to record, skip it
+                const string particleName = track.GetParticleName().Data();
+                if (find(particlesToRecord.begin(), particlesToRecord.end(), particleName) ==
+                    particlesToRecord.end()) {
+                    continue;
+                }
+                // record energy at the point of exit from the "detector" sphere (if it exits)
+                const auto& hits = track.GetHits();
+                for (int i = 0; i < hits.GetNumberOfHits(); i++) {
+                    if (hits.GetVolumeName(i) == detectorName && hits.GetProcessName(i) == "Transportation") {
+                        // We can use either the kinetic energy at the point of exit from the detector or the
+                        // initial energy depending on our goals energies[particleName] =
+                        // hits.GetKineticEnergy(i);
+                        energies[particleName] = track.GetInitialKineticEnergy();
+                        trees[particleName]->Fill();
+                    }
+                }
+            }
+        }
+    }
+
+    outputFile->Write();
+    outputFile->Close();
+
+    /*
+     * Histograms can be merged into a single file using the following root code:
+
+TFile* U238 = TFile::Open("U238/gammas.root");
+TFile* Th232 = TFile::Open("Th232/gammas.root");
+
+TH1D* U238h = U238->Get<TH1D>("histogram");
+TH1D* Th232h = Th232->Get<TH1D>("histogram");
+
+// rename histograms
+U238h->SetName("U238");
+Th232h->SetName("Th232");
+
+// write to file
+TFile* outputFile = new TFile("merged.root", "RECREATE");
+
+U238h->Write();
+Th232h->Write();
+
+outputFile->Close();
+
+     */
+}

examples/15.DecaySpectra/TreeToHistogram.C

@@ -0,0 +1,23 @@
+
+using namespace std;
+
+void TreeToHistogram(const char* filename) {
+    TFile* file = TFile::Open(filename);
+
+    TTree* tree = nullptr;
+    tree = file->Get<TTree>("gammas");
+
+    TH1D* histogram = new TH1D("histogram", "histogram", 10000, 0, 5000);
+
+    double energy;
+    tree->SetBranchAddress("energy", &energy);
+    for (int i = 0; i < tree->GetEntries(); i++) {
+        tree->GetEntry(i);
+        histogram->Fill(energy);
+    }
+    // write to file
+    TFile* outputFile = new TFile("gammas.root", "RECREATE");
+
+    histogram->Write();
+    outputFile->Close();
+}

examples/15.DecaySpectra/geometry/setup.gdml

@@ -0,0 +1,44 @@
+<?xml version="1.0" encoding="utf-8" standalone="no" ?>
+
+<!DOCTYPE gdml [
+        <!ENTITY materials SYSTEM "https://rest-for-physics.github.io/materials/materials.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>
+        <constant name="worldSize" value="10"/>
+    </define>
+
+    &materials;
+
+    <solids>
+        <box name="worldSolid" x="worldSize" y="worldSize" z="worldSize" lunit="mm"/>
+        <orb name="detectorSphereSolid" r="2.5" lunit="mm"/>
+    </solids>
+
+    <structure>
+
+        <volume name="detectorVolume">
+            <materialref ref="G4_Pb"/>
+            <solidref ref="detectorSphereSolid"/>
+        </volume>
+
+        <volume name="World">
+            <materialref ref="Vacuum"/>
+            <solidref ref="worldSolid"/>
+
+            <physvol name="detector">
+                <volumeref ref="detectorVolume"/>
+                <position name="detectorPosition" unit="mm" x="0" y="0" z="0"/>
+            </physvol>
+        </volume>
+
+    </structure>
+
+    <setup name="Default" version="1.0">
+        <world ref="World"/>
+    </setup>
+
+</gdml>

examples/15.DecaySpectra/simulation.rml

@@ -0,0 +1,83 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no" ?>
+
+<restG4>
+
+    <globals>
+        <parameter name="mainDataPath" value=""/>
+        <parameter name="verboseLevel" value="essential"/>
+        <variable name="REST_ISOTOPE" value="U238"/>
+    </globals>
+
+    <TRestRun name="Metadata">
+        <parameter name="experimentName" value="DecaySpectra"/>
+        <parameter name="runType" value="simulation"/>
+        <parameter name="runNumber" value="1"/>
+        <parameter name="runTag" value="${REST_ISOTOPE}"/>
+        <parameter name="outputFileName" value="Run[fRunNumber]_[fRunTag]_[fExperimentName].root"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="overwrite" value="off"/>
+        <parameter name="readOnly" value="false"/>
+    </TRestRun>
+
+    <TRestGeant4Metadata name="DecaySpectra" title="DecaySpectra">
+
+        <parameter name="gdmlFile" value="geometry/setup.gdml"/>
+        <parameter name="subEventTimeDelay" value="100" units="us"/>
+
+        <parameter name="seed" value="1000"/>
+
+        <!-- The number of events to be simulated is now defined in TRestGeant4Metadata -->
+        <parameter name="nEvents" value="10000000"/>
+        <parameter name="saveAllEvents" value="true"/>
+
+        <generator type="point" position="(0,0,0)" units="mm">
+            <source particle="${REST_ISOTOPE}" fullChain="on">
+                <angular type="isotropic"/>
+                <energy type="mono" energy="0keV"/>
+            </source>
+        </generator>
+
+        <detector>
+            <parameter name="energyRange" value="(0,1)" units="GeV"/>
+            <volume name="detector" sensitive="true" maxStepSize="1mm"/>
+        </detector>
+
+    </TRestGeant4Metadata>
+
+    <TRestGeant4PhysicsLists name="default" title="First physics list implementation.">
+
+        <parameter name="ionLimitStepList" value="F20,Ne20"/>
+
+        <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"/>
+
+        <!-- EM Physics lists -->
+        <physicsList name="G4EmStandardPhysics_option4"> <!-- "G4EmPenelopePhysics", "G4EmStandardPhysics_option3" -->
+            <option name="pixe" value="true"/>
+        </physicsList>
+
+        <!-- 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>

src/PhysicsList.cxx

@@ -335,7 +335,7 @@ void PhysicsList::ConstructProcess() {
                 if (fRestPhysicsLists->GetIonStepList()[n] == G4IonTable::GetIonTable()->GetIonName(Z, A)) {
                     G4ParticleDefinition* particle = G4IonTable::GetIonTable()->GetIon(Z, A, 0);
                     G4String particle_name = G4IonTable::GetIonTable()->GetIonName(Z, A, 0);
-                    cout << "Found ion: " << particle_name << " Z " << Z << " A " << A << endl;
+                    RESTDebug << "Found ion: " << particle_name << " Z " << Z << " A " << A << RESTendl;
                     G4ProcessManager* processManager = particle->GetProcessManager();
                     processManager->AddDiscreteProcess(new G4StepLimiter("ionStep"));
                 }