From 0a63b4ed69fcdc9631d7f2cfb2e31cfd97a5f689 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Thu, 13 Apr 2023 12:54:56 +0200 Subject: [PATCH 1/8] add new example to use restG4 to generate spectra of decay particles as a root TTree --- examples/15.DecaySpectra/GenerateTrees.C | 64 +++++++++++++++ examples/15.DecaySpectra/geometry/setup.gdml | 44 +++++++++++ examples/15.DecaySpectra/simulation.rml | 83 ++++++++++++++++++++ 3 files changed, 191 insertions(+) create mode 100644 examples/15.DecaySpectra/GenerateTrees.C create mode 100644 examples/15.DecaySpectra/geometry/setup.gdml create mode 100644 examples/15.DecaySpectra/simulation.rml diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTrees.C new file mode 100644 index 00000000..a91ce40d --- /dev/null +++ b/examples/15.DecaySpectra/GenerateTrees.C @@ -0,0 +1,64 @@ + +using namespace std; + +void GenerateTrees(const char *filename) { + gSystem->Load("libRestFramework.so"); + gSystem->Load("libRestGeant4.so"); + + TRestRun run(filename); + cout << "Number of entries: " << run.GetEntries() << endl; + + TRestGeant4Event *event = run.GetInputEvent(); + run.GetEntry(0); + + const TString primaryName = event->GetPrimaryEventParticleName(); + // 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 particlesToRecord = {"gamma", "e-", "e+", "alpha"}; + // 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())); + } + + // create a branch for each particle + map energies; + for (const auto &particle: particlesToRecord) { + energies[particle] = 0; + trees[particle]->Branch("energy", &energies[particle], "energy/D"); + } + + const string detectorName = "detector"; + 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") { + // this is more accurate than track.GetInitialKineticEnergy() + // because it takes into account energy loss in the detector + energies[particleName] = hits.GetKineticEnergy(i); + trees[particleName]->Fill(); + } + } + } + } + + f->Write(); + f->Close(); +} diff --git a/examples/15.DecaySpectra/geometry/setup.gdml b/examples/15.DecaySpectra/geometry/setup.gdml new file mode 100644 index 00000000..b03e37b9 --- /dev/null +++ b/examples/15.DecaySpectra/geometry/setup.gdml @@ -0,0 +1,44 @@ + + + + ]> + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/15.DecaySpectra/simulation.rml b/examples/15.DecaySpectra/simulation.rml new file mode 100644 index 00000000..060eb9a1 --- /dev/null +++ b/examples/15.DecaySpectra/simulation.rml @@ -0,0 +1,83 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 708cef0d8003a1befc941fa500e26c7e38102115 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 13 Apr 2023 11:00:30 +0000 Subject: [PATCH 2/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- examples/15.DecaySpectra/GenerateTrees.C | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTrees.C index a91ce40d..b8e90343 100644 --- a/examples/15.DecaySpectra/GenerateTrees.C +++ b/examples/15.DecaySpectra/GenerateTrees.C @@ -1,24 +1,24 @@ using namespace std; -void GenerateTrees(const char *filename) { +void GenerateTrees(const char* filename) { gSystem->Load("libRestFramework.so"); gSystem->Load("libRestGeant4.so"); TRestRun run(filename); cout << "Number of entries: " << run.GetEntries() << endl; - TRestGeant4Event *event = run.GetInputEvent(); + TRestGeant4Event* event = run.GetInputEvent(); run.GetEntry(0); const TString primaryName = event->GetPrimaryEventParticleName(); // Create a new tree to store gammas from the primary particle - TFile *f = new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE"); + TFile* f = new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE"); - const vector particlesToRecord = {"gamma", "e-", "e+", "alpha"}; + const vector particlesToRecord = {"gamma", "e-", "e+", "alpha"}; // create a new tree for each particle - map < string, TTree * > trees; - for (const auto &particle: particlesToRecord) { + map 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())); @@ -26,7 +26,7 @@ void GenerateTrees(const char *filename) { // create a branch for each particle map energies; - for (const auto &particle: particlesToRecord) { + for (const auto& particle : particlesToRecord) { energies[particle] = 0; trees[particle]->Branch("energy", &energies[particle], "energy/D"); } @@ -39,7 +39,7 @@ void GenerateTrees(const char *filename) { } run.GetEntry(n); // iterate over all particles in the event - for (const auto &track: event->GetTracks()) { + 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) == @@ -47,7 +47,7 @@ void GenerateTrees(const char *filename) { continue; } // record energy at the point of exit from the "detector" sphere (if it exits) - const auto &hits = track.GetHits(); + const auto& hits = track.GetHits(); for (int i = 0; i < hits.GetNumberOfHits(); i++) { if (hits.GetVolumeName(i) == detectorName && hits.GetProcessName(i) == "Transportation") { // this is more accurate than track.GetInitialKineticEnergy() From 30d20325ce90b3d8c4a24b4446b90b7d8831f6a7 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Tue, 25 Apr 2023 12:23:24 +0200 Subject: [PATCH 3/8] move print to debug --- src/PhysicsList.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PhysicsList.cxx b/src/PhysicsList.cxx index 55748595..23c4d4cf 100644 --- a/src/PhysicsList.cxx +++ b/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")); } From 3c138c304dcd1fc2e40394fe1c1f620bf9481347 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Tue, 2 May 2023 10:09:20 +0200 Subject: [PATCH 4/8] use solid detector --- examples/15.DecaySpectra/geometry/setup.gdml | 4 ++-- examples/15.DecaySpectra/simulation.rml | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/15.DecaySpectra/geometry/setup.gdml b/examples/15.DecaySpectra/geometry/setup.gdml index b03e37b9..0c10079a 100644 --- a/examples/15.DecaySpectra/geometry/setup.gdml +++ b/examples/15.DecaySpectra/geometry/setup.gdml @@ -15,13 +15,13 @@ - + - + diff --git a/examples/15.DecaySpectra/simulation.rml b/examples/15.DecaySpectra/simulation.rml index 060eb9a1..441bb6d3 100644 --- a/examples/15.DecaySpectra/simulation.rml +++ b/examples/15.DecaySpectra/simulation.rml @@ -64,19 +64,19 @@ - - + + - - - - - + + + + + From e4566d29dd42c4158df94aa94ca5e9606f19bb1c Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 3 May 2023 09:19:57 +0200 Subject: [PATCH 5/8] update macro and settings --- examples/15.DecaySpectra/GenerateTrees.C | 18 +++++++++++++++--- examples/15.DecaySpectra/geometry/setup.gdml | 4 ++-- examples/15.DecaySpectra/simulation.rml | 2 +- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTrees.C index b8e90343..837b7cf4 100644 --- a/examples/15.DecaySpectra/GenerateTrees.C +++ b/examples/15.DecaySpectra/GenerateTrees.C @@ -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 particlesToRecord = {"gamma", "e-", "e+", "alpha"}; + // store particles and name of the root TTree + const map particleToString = { + {"gamma", "photons"}, {"e-", "electrons"}, {"e+", "positrons"}, + {"proton", "protons"}, {"neutron", "neutrons"}, {"alpha", "alphas"}, + }; + vector particlesToRecord; + for (const auto& particle : particleToString) { + particlesToRecord.push_back(particle.first); + } // create a new tree for each particle map 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 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) { diff --git a/examples/15.DecaySpectra/geometry/setup.gdml b/examples/15.DecaySpectra/geometry/setup.gdml index 0c10079a..1bb0aec0 100644 --- a/examples/15.DecaySpectra/geometry/setup.gdml +++ b/examples/15.DecaySpectra/geometry/setup.gdml @@ -8,14 +8,14 @@ xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd"> - + &materials; - + diff --git a/examples/15.DecaySpectra/simulation.rml b/examples/15.DecaySpectra/simulation.rml index 441bb6d3..2621f08d 100644 --- a/examples/15.DecaySpectra/simulation.rml +++ b/examples/15.DecaySpectra/simulation.rml @@ -20,7 +20,7 @@ - + From fb26ec606c1df30ad04e70eaac1ae61eb9d49d5e Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 3 May 2023 12:56:07 +0200 Subject: [PATCH 6/8] support multiple files --- examples/15.DecaySpectra/GenerateTrees.C | 141 +++++++++++++---------- 1 file changed, 83 insertions(+), 58 deletions(-) diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTrees.C index 837b7cf4..06f44bcd 100644 --- a/examples/15.DecaySpectra/GenerateTrees.C +++ b/examples/15.DecaySpectra/GenerateTrees.C @@ -1,76 +1,101 @@ +#include +#include using namespace std; -void GenerateTrees(const char* filename) { +void GenerateTrees(const char *filenameOrDirectoryWithFiles) { gSystem->Load("libRestFramework.so"); gSystem->Load("libRestGeant4.so"); - TRestRun run(filename); - cout << "Number of entries: " << run.GetEntries() << endl; + vector filenames; + if (filesystem::is_directory(filenameOrDirectoryWithFiles)) { + for (const auto &entry: filesystem::directory_iterator(filenameOrDirectoryWithFiles)) { + filenames.push_back(entry.path()); + } + } else { + filenames.push_back(filenameOrDirectoryWithFiles); + } - TRestGeant4Event* event = run.GetInputEvent(); - run.GetEntry(0); + bool first = true; + for (const auto &filename: filenames) { + cout << "Processing file: " << filename << endl; - const TString primaryName = event->GetPrimaryEventParticleName(); - // Create a new tree to store gammas from the primary particle - TFile* f = new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE"); + TRestRun run(filename); + cout << "Number of entries: " << run.GetEntries() << endl; - // store particles and name of the root TTree - const map particleToString = { - {"gamma", "photons"}, {"e-", "electrons"}, {"e+", "positrons"}, - {"proton", "protons"}, {"neutron", "neutrons"}, {"alpha", "alphas"}, - }; - vector particlesToRecord; - for (const auto& particle : particleToString) { - particlesToRecord.push_back(particle.first); - } - // create a new tree for each particle - map 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 - 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()); - } + TRestGeant4Event *event = run.GetInputEvent(); + run.GetEntry(0); - // create a branch for each particle - map energies; - for (const auto& particle : particlesToRecord) { - cout << "Creating branch for " << particle << endl; - energies[particle] = 0; - trees[particle]->Branch("energy", &energies[particle], "energy/D"); - } + TFile *outputFile = nullptr; + const map particleToString = { + {"gamma", "gammas"}, + /* + {"e-", "electrons"}, + {"e+", "positrons"}, + {"proton", "protons"}, + {"neutron", "neutrons"}, + {"alpha", "alphas"}, + */ + }; + vector particlesToRecord; + for (const auto &particle: particleToString) { + particlesToRecord.push_back(particle.first); + } + map > trees; + map energies; + + 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"); - 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; + // 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] = make_unique(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; } - 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; + + 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; } - // 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") { - // this is more accurate than track.GetInitialKineticEnergy() - // because it takes into account energy loss in the detector - energies[particleName] = hits.GetKineticEnergy(i); - trees[particleName]->Fill(); + 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(); + } } } } - } - f->Write(); - f->Close(); + outputFile->Write(); + outputFile->Close(); + } } From 07b99010569fda3a215af764f22319803cb8b762 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 3 May 2023 10:56:26 +0000 Subject: [PATCH 7/8] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- examples/15.DecaySpectra/GenerateTrees.C | 53 ++++++++++++------------ 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTrees.C index 06f44bcd..b30f5597 100644 --- a/examples/15.DecaySpectra/GenerateTrees.C +++ b/examples/15.DecaySpectra/GenerateTrees.C @@ -3,13 +3,13 @@ using namespace std; -void GenerateTrees(const char *filenameOrDirectoryWithFiles) { +void GenerateTrees(const char* filenameOrDirectoryWithFiles) { gSystem->Load("libRestFramework.so"); gSystem->Load("libRestGeant4.so"); - vector filenames; + vector filenames; if (filesystem::is_directory(filenameOrDirectoryWithFiles)) { - for (const auto &entry: filesystem::directory_iterator(filenameOrDirectoryWithFiles)) { + for (const auto& entry : filesystem::directory_iterator(filenameOrDirectoryWithFiles)) { filenames.push_back(entry.path()); } } else { @@ -17,47 +17,48 @@ void GenerateTrees(const char *filenameOrDirectoryWithFiles) { } bool first = true; - for (const auto &filename: filenames) { + 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* event = run.GetInputEvent(); run.GetEntry(0); - TFile *outputFile = nullptr; - const map particleToString = { - {"gamma", "gammas"}, - /* - {"e-", "electrons"}, - {"e+", "positrons"}, - {"proton", "protons"}, - {"neutron", "neutrons"}, - {"alpha", "alphas"}, - */ + TFile* outputFile = nullptr; + const map particleToString = { + {"gamma", "gammas"}, + /* + {"e-", "electrons"}, + {"e+", "positrons"}, + {"proton", "protons"}, + {"neutron", "neutrons"}, + {"alpha", "alphas"}, + */ }; - vector particlesToRecord; - for (const auto &particle: particleToString) { + vector particlesToRecord; + for (const auto& particle : particleToString) { particlesToRecord.push_back(particle.first); } - map > trees; + map> trees; map energies; 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"); + outputFile = + new TFile(TString::Format("ParticlesFrom%sDecay.root", primaryName.Data()), "RECREATE"); // create a new tree for each particle - for (const auto &particle: particlesToRecord) { + 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] = make_unique(particleToString.at(particle).c_str(), particleToString.at(particle).c_str()); } // create a branch for each particle - for (const auto &particle: particlesToRecord) { + for (const auto& particle : particlesToRecord) { energies[particle] = 0; trees[particle]->Branch("energy", &energies[particle], "energy/D"); } @@ -74,7 +75,7 @@ void GenerateTrees(const char *filenameOrDirectoryWithFiles) { } run.GetEntry(n); // iterate over all particles in the event - for (const auto &track: event->GetTracks()) { + 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) == @@ -82,12 +83,12 @@ void GenerateTrees(const char *filenameOrDirectoryWithFiles) { continue; } // record energy at the point of exit from the "detector" sphere (if it exits) - const auto &hits = track.GetHits(); + 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); + // 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(); } From fd0a19fc5fcca940f4e60cda74bb334bc4755db7 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Thu, 4 May 2023 10:54:10 +0200 Subject: [PATCH 8/8] update tree generation code, add code to produce histogram --- .../{GenerateTrees.C => GenerateTree.C} | 72 ++++++++++++------- examples/15.DecaySpectra/TreeToHistogram.C | 23 ++++++ 2 files changed, 71 insertions(+), 24 deletions(-) rename examples/15.DecaySpectra/{GenerateTrees.C => GenerateTree.C} (71%) create mode 100644 examples/15.DecaySpectra/TreeToHistogram.C diff --git a/examples/15.DecaySpectra/GenerateTrees.C b/examples/15.DecaySpectra/GenerateTree.C similarity index 71% rename from examples/15.DecaySpectra/GenerateTrees.C rename to examples/15.DecaySpectra/GenerateTree.C index b30f5597..a2f20ccb 100644 --- a/examples/15.DecaySpectra/GenerateTrees.C +++ b/examples/15.DecaySpectra/GenerateTree.C @@ -3,7 +3,7 @@ using namespace std; -void GenerateTrees(const char* filenameOrDirectoryWithFiles) { +void GenerateTree(const char* filenameOrDirectoryWithFiles) { gSystem->Load("libRestFramework.so"); gSystem->Load("libRestGeant4.so"); @@ -11,11 +11,30 @@ void GenerateTrees(const char* filenameOrDirectoryWithFiles) { 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 particleToString = { + {"gamma", "gammas"}, + /* + {"e-", "electrons"}, + {"e+", "positrons"}, + {"proton", "protons"}, + {"neutron", "neutrons"}, + {"alpha", "alphas"}, + */ + }; + vector particlesToRecord; + for (const auto& particle : particleToString) { + particlesToRecord.push_back(particle.first); + } + map trees; + map energies; + bool first = true; for (const auto& filename : filenames) { cout << "Processing file: " << filename << endl; @@ -26,24 +45,6 @@ void GenerateTrees(const char* filenameOrDirectoryWithFiles) { TRestGeant4Event* event = run.GetInputEvent(); run.GetEntry(0); - TFile* outputFile = nullptr; - const map particleToString = { - {"gamma", "gammas"}, - /* - {"e-", "electrons"}, - {"e+", "positrons"}, - {"proton", "protons"}, - {"neutron", "neutrons"}, - {"alpha", "alphas"}, - */ - }; - vector particlesToRecord; - for (const auto& particle : particleToString) { - particlesToRecord.push_back(particle.first); - } - map> trees; - map energies; - if (first) { const TString primaryName = event->GetPrimaryEventParticleName(); // Create a new tree to store gammas from the primary particle @@ -53,8 +54,8 @@ void GenerateTrees(const char* filenameOrDirectoryWithFiles) { // 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] = make_unique(particleToString.at(particle).c_str(), - particleToString.at(particle).c_str()); + trees[particle] = + new TTree(particleToString.at(particle).c_str(), particleToString.at(particle).c_str()); } // create a branch for each particle @@ -95,8 +96,31 @@ void GenerateTrees(const char* filenameOrDirectoryWithFiles) { } } } - - outputFile->Write(); - outputFile->Close(); } + + 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("histogram"); +TH1D* Th232h = Th232->Get("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(); + + */ } diff --git a/examples/15.DecaySpectra/TreeToHistogram.C b/examples/15.DecaySpectra/TreeToHistogram.C new file mode 100644 index 00000000..3754c4f7 --- /dev/null +++ b/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("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(); +}