Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New example to generate decay spectra as a tree #103

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
Prev Previous commit
Next Next commit
update macro and settings
  • Loading branch information
lobis committed May 3, 2023
commit e4566d29dd42c4158df94aa94ca5e9606f19bb1c
18 changes: 15 additions & 3 deletions examples/15.DecaySpectra/GenerateTrees.C
Original file line number Diff line number Diff line change
@@ -15,23 +15,35 @@ void GenerateTrees(const char* filename) {
// Create a new tree to store gammas from the primary particle
TFile* f = new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE");

const vector<string> particlesToRecord = {"gamma", "e-", "e+", "alpha"};
// store particles and name of the root TTree
const map<string, string> particleToString = {
{"gamma", "photons"}, {"e-", "electrons"}, {"e+", "positrons"},
{"proton", "protons"}, {"neutron", "neutrons"}, {"alpha", "alphas"},
};
vector<string> particlesToRecord;
for (const auto& particle : particleToString) {
particlesToRecord.push_back(particle.first);
}
// create a new tree for each particle
map<string, TTree*> trees;
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(TString::Format("decay_%s", particle.c_str()),
TString::Format("%s particles", particle.c_str()));
cout << "Creating tree for " << particle << " with name " << particleToString.at(particle) << endl;
trees[particle] =
new TTree(particleToString.at(particle).c_str(), particleToString.at(particle).c_str());
}

// create a branch for each particle
map<string, double> energies;
for (const auto& particle : particlesToRecord) {
cout << "Creating branch for " << particle << endl;
energies[particle] = 0;
trees[particle]->Branch("energy", &energies[particle], "energy/D");
}

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) {
4 changes: 2 additions & 2 deletions examples/15.DecaySpectra/geometry/setup.gdml
Original file line number Diff line number Diff line change
@@ -8,14 +8,14 @@
xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">

<define>
<constant name="worldSize" value="100"/>
<constant name="worldSize" value="10"/>
</define>

&materials;

<solids>
<box name="worldSolid" x="worldSize" y="worldSize" z="worldSize" lunit="mm"/>
<orb name="detectorSphereSolid" r="5" lunit="mm"/>
<orb name="detectorSphereSolid" r="2.5" lunit="mm"/>
</solids>

<structure>
2 changes: 1 addition & 1 deletion examples/15.DecaySpectra/simulation.rml
Original file line number Diff line number Diff line change
@@ -20,7 +20,7 @@
<parameter name="readOnly" value="false"/>
</TRestRun>

<TRestGeant4Metadata name="DecaySpectra" title="DecaySpectra_{REST_ISOTOPE}">
<TRestGeant4Metadata name="DecaySpectra" title="DecaySpectra">

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