BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
HistogramDefSet.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "HistogramDef.hh"
20#include "HistogramDefSet.hh"
21#include "RBDSException.hh"
22#include "SpectraParticles.hh"
23
24#include <map>
25#include <regex>
26#include <set>
27#include <stdexcept>
28#include <string>
29
30HistogramDefSet::HistogramDefSet(const std::string& branchNameIn,
31 const HistogramDef* baseDefinitionIn,
32 const std::set<ParticleSpec>& particlesSpecs,
33 const std::string& particleSpecificationIn):
34 branchName(branchNameIn),
35 dynamicallyStoreIons(false),
36 dynamicallyStoreParticles(particlesSpecs.empty()),
37 what(writewhat::all),
38 topN(1)
39{
40 if (!baseDefinitionIn)
41 {throw std::invalid_argument("invalid histogram definition");}
42
43 baseDefinition = baseDefinitionIn->Clone();
44
45 if (!particlesSpecs.empty())
46 {
47 for (const auto& particleSpec : particlesSpecs)
48 {
49 HistogramDef* h = baseDefinitionIn->Clone();
50 std::string thisHistName = "Spectra_" + h->histName + "_";
51 if (particleSpec.first > 0) // prepend + character for positive values so we always have + or -
52 {thisHistName += "+";}
53 thisHistName += std::to_string(particleSpec.first);
54 h->histName = thisHistName;
55 std::string suffix;
56 switch (particleSpec.second)
57 {
58 case RBDS::SpectraParticles::primary:
59 {suffix = "_Primary"; break;}
60 case RBDS::SpectraParticles::secondary:
61 {suffix = "_Secondary"; break;}
62 default:
63 {break;}
64 }
65 h->histName += suffix;
66 h->selection = AddPDGFilterToSelection(particleSpec, h->selection, branchName);
67 definitions[particleSpec] = h;
68 definitionsV.push_back(h);
69 }
70 }
71 else
72 {
73 std::string spec = particleSpecificationIn;
74 std::transform(spec.begin(), spec.end(), spec.begin(), ::tolower);
75 spec = RemoveSubString(spec, "{");
76 spec = RemoveSubString(spec, "}");
77
78 if (spec.find("all") != std::string::npos || spec.find("ions") != std::string::npos)
79 {dynamicallyStoreIons = true;}
80 if (spec.find("all") != std::string::npos || spec.find("particles") != std::string::npos)
81 {dynamicallyStoreParticles = true;}
82
83 std::map<std::string, writewhat> keys = {{"all", writewhat::all},
84 {"particles", writewhat::particles},
85 {"ions", writewhat::ions}};
86
87 auto search = keys.find(spec);
88 if (search != keys.end())
89 {what = search->second;}
90 else
91 {
92 std::regex topNR("top(\\d+)(?:particles|ions)*");
93 std::smatch match;
94 if (std::regex_search(spec, match, topNR))
95 {
96 try
97 {topN = std::stoi(match[1]);}
98 catch (...)
99 {topN = 1;}
100 if (spec.find("particles") != std::string::npos)
101 {what = writewhat::topNParticles;}
102 else if (spec.find("ions") != std::string::npos)
103 {what = writewhat::topNIons;}
104 else
105 {what = writewhat::topN;}
106 }
107 else // this will happen if 'top' isn't in the specification or it generally doesn't match
108 {throw RBDSException("Invalid particle specifier \"" + particleSpecificationIn + "\"");}
109 }
110 }
111}
112
113HistogramDefSet::~HistogramDefSet()
114{
115 delete baseDefinition;
116 for (auto kv : definitions)
117 {delete kv.second;}
118}
119
120std::string HistogramDefSet::RemoveSubString(const std::string& stringIn,
121 const std::string& wordToRemove) const
122{
123 std::string result = stringIn;
124 while (result.find(wordToRemove) != std::string::npos)
125 {
126 size_t pos = result.find(wordToRemove);
127 result.erase(pos, wordToRemove.size());
128 }
129 return result;
130}
131
132std::ostream& operator<< (std::ostream &out, const HistogramDefSet& s)
133{
134 out << "Spectra: " << s.baseDefinition->histName << "\n";
135 for (const auto* d : s.definitionsV)
136 {out << *d;}
137 return out;
138}
139
140std::string HistogramDefSet::AddPDGFilterToSelection(const ParticleSpec& particleSpec,
141 const std::string& selection,
142 const std::string& branchName)
143{
144 long long int pdgID = particleSpec.first;
145
146 if (pdgID == 0) // '0' is our code for total or really all particles in 1x histogram
147 {return selection;} // just return whatever selection is
148
149 RBDS::SpectraParticles flag = particleSpec.second;
150 std::string flagFilter;
151 switch (flag)
152 {
153 case RBDS::SpectraParticles::primary:
154 {flagFilter = "&&" + branchName + ".parentID==0"; break;}
155 case RBDS::SpectraParticles::secondary:
156 {flagFilter = "&&" + branchName + ".parentID>0"; break;}
157 default:
158 {break;}
159 }
160 std::string filter = branchName+".partID=="+std::to_string(pdgID) + flagFilter;
161
162 // input selection could be:
163 // 1
164 // samplerName.weight
165 // samplerName.variable*number/number
166 // samplerName.variable<value // Boolean expression
167 // variable*(Boolean)
168 // we need to put in the bonus Boolean if needed
169 // check if it has a Boolean expression already in it
170 std::string result;
171 std::regex boolOperator("&&|[<>!=]=|[<>]|\\|\\|");
172 std::smatch match;
173 if (std::regex_search(selection, match, boolOperator))
174 {// has boolean operator somewhere in it
175 std::string afterBool = match.suffix();
176 std::size_t bracketPos = afterBool.find(')');
177 result = selection; // copy it
178 if (bracketPos == std::string::npos) // no bracket found - so just Boolean condition on its own
179 {result += "&&" + filter;}
180 else
181 {result.insert(match.position()+ match.length() + bracketPos, "&&"+filter);}
182 }
183 else if (selection.empty() || selection == "1") // technically the selection shouldn't be empty, but just in case...
184 {result = filter;}
185 else
186 {result = selection + "*("+filter+")";}
187 return result;
188}
189
Specification for a set of histograms.
std::vector< HistogramDef * > definitionsV
Vector version for easy iteration.
std::string RemoveSubString(const std::string &stringIn, const std::string &wordToRemove) const
Remove a substring from a string.
Common specification for a histogram.
Definition: HistogramDef.hh:34
virtual HistogramDef * Clone() const =0
Copy this instance. Virtual to be overridden in derived classes.
General exception with possible name of object and message.