19#include "BDSArray1DCoords.hh"
20#include "BDSArray2DCoords.hh"
21#include "BDSArray3DCoords.hh"
22#include "BDSArray4DCoords.hh"
24#include "BDSException.hh"
25#include "BDSFieldLoaderBDSIM.hh"
26#include "BDSFieldValue.hh"
27#include "BDSUtilities.hh"
32#include "CLHEP/Units/SystemOfUnits.h"
46#include "src-external/gzstream/gzstream.h"
54 headerMustBePositiveKeys({
"nx",
"ny",
"nz",
"nt"}),
55 indexOfFirstFieldValue(0)
58 {BDSDimensionType::x, {
"nx",
"xmin",
"xmax"}},
59 {BDSDimensionType::y, {
"ny",
"ymin",
"ymax"}},
60 {BDSDimensionType::z, {
"nz",
"zmin",
"zmax"}},
61 {BDSDimensionType::t, {
"nt",
"tmin",
"tmax"}}
75 std::vector<G4String> allKeys = {
"nx",
"ny",
"nz",
"nt",
76 "xmin",
"xmax",
"ymin",
"ymax",
77 "zmin",
"zmax",
"tmin",
"tmax"};
78 for (
auto& key : allKeys)
80 for (
auto& key : headerMustBePositiveKeys)
84 indexOfFirstFieldValue = 0;
124 const unsigned int nDim)
126 G4String functionName =
"BDSIM Field Format> ";
130 long long unsigned int currentLineNumber = 0;
133 bool validFile = file.rdbuf()->is_open();
135 bool validFile = file.is_open();
138 {
throw BDSException(__METHOD_NAME__,
"Invalid file name or no such file named \"" + fileName +
"\"");}
140 {G4cout << functionName <<
"Loading \"" << fileName <<
"\"" << G4endl;}
143 unsigned long xIndex = 0;
144 unsigned long yIndex = 0;
145 unsigned long zIndex = 0;
146 G4bool intoData =
false;
152 G4int n1 = 1, n2 = 1, n3 = 1, n4 = 1;
155 G4String nominalOrder =
"xyzt";
157 float maximumFieldValue = 0;
158 float minimumFieldValue = 0;
162 std::getline(file, line);
166 if (std::all_of(line.begin(), line.end(), isspace))
169 if (currentLineNumber > 100)
171 G4String msg = functionName+
"did not encounter column specification line starting with '!' in the first 100 lines\n";
172 msg +=
"This does not appear to be a BDSIM-format field map.";
186 std::smatch matchHeaderNumber;
187 std::regex keyValue(R
"((\w*)\s*>\s*([0-9eE.+-]+))");
188 if (std::regex_search(line, matchHeaderNumber, keyValue))
190 if (matchHeaderNumber.size() < 2)
191 {Terminate(functionName+
"Invalid key definition in field format: " + line);}
194 G4String key = G4String(matchHeaderNumber[1]);
198 if (header.find(key) == header.end())
199 {Terminate(functionName+
"Invalid key \"" + key +
"\" in header");}
203 {value = std::stod(matchHeaderNumber[2]);}
204 catch (
const std::invalid_argument&)
205 {Terminate(functionName+
"Invalid argument " + std::string(matchHeaderNumber[2]));}
206 catch (
const std::out_of_range&)
207 {Terminate(functionName+
"Number out of range " + std::string(matchHeaderNumber[2]));}
209 if (std::find(headerMustBePositiveKeys.begin(), headerMustBePositiveKeys.end(), key) != headerMustBePositiveKeys.end())
212 {Terminate(functionName+
"Number of points in dimension must be greater than 0 -> see \"" + key +
"\"");}
220 std::smatch matchHeaderString;
224 std::regex keyWord(R
"((\w+)\s*>\s*([a-zA-Z]{1,4})\b\s*)");
225 if (std::regex_search(line, matchHeaderString, keyWord))
227 loopOrder = G4String(matchHeaderString[2]);
229 G4bool test1 = loopOrder ==
"xyzt";
230 G4bool test2 = loopOrder ==
"tzyx";
231 if (! (test1 || test2) )
232 {Terminate(functionName+
"loopOrder> header variable is not one of (\"xyzt\" or \"tzyx\")");}
239 std::regex columnRow(
"^\\s*!");
240 if (std::regex_search(line, columnRow))
244 std::regex afterExclamation(
"\\s*!\\s*(.+)");
246 std::regex_search(line, match, afterExclamation);
247 std::string restOfLine = match[1];
248 std::string columnName;
249 std::istringstream restOfLineSS(restOfLine);
250 std::vector<G4String> columnNames;
251 while (restOfLineSS >> columnName)
254 if (columnName.find(
"Fx") != std::string::npos)
255 {xIndex = nColumns;
continue;}
256 else if (columnName.find(
"Fy") != std::string::npos)
257 {yIndex = nColumns;
continue;}
258 else if (columnName.find(
"Fz") != std::string::npos)
259 {zIndex = nColumns;
continue;}
261 {columnNames.emplace_back(columnName);}
263 lineData.resize(nColumns + 1);
264 indexOfFirstFieldValue = std::min({xIndex, yIndex, zIndex});
266 for (
const auto& key : headerMustBePositiveKeys)
269 {Terminate(functionName+
"Number of points in dimension must be greater than 0 -> see \"" + key +
"\"");}
272 if (nColumns < (nDim + 3))
273 {Terminate(functionName+
"Too few columns for " + std::to_string(nDim) +
"D field loading");}
281 auto keys = dimKeyMap[firstDim];
282 n1 = G4int(header[keys.number]);
284 header[keys.min] * CLHEP::cm,
285 header[keys.max] * CLHEP::cm,
293 auto fKeys = dimKeyMap[firstDim];
294 auto sKeys = dimKeyMap[secondDim];
295 n1 = G4int(header[fKeys.number]);
296 n2 = G4int(header[sKeys.number]);
298 header[fKeys.min] * CLHEP::cm,
299 header[fKeys.max] * CLHEP::cm,
300 header[sKeys.min] * CLHEP::cm,
301 header[sKeys.max] * CLHEP::cm,
311 auto fKeys = dimKeyMap[firstDim];
312 auto sKeys = dimKeyMap[secondDim];
313 auto tKeys = dimKeyMap[thirdDim];
314 n1 = G4int(header[fKeys.number]);
315 n2 = G4int(header[sKeys.number]);
316 n3 = G4int(header[tKeys.number]);
318 header[fKeys.min] * CLHEP::cm,
319 header[fKeys.max] * CLHEP::cm,
320 header[sKeys.min] * CLHEP::cm,
321 header[sKeys.max] * CLHEP::cm,
322 header[tKeys.min] * CLHEP::cm,
323 header[tKeys.max] * CLHEP::cm,
331 n1 = G4int(header[
"nx"]);
332 n2 = G4int(header[
"ny"]);
333 n3 = G4int(header[
"nz"]);
334 n4 = G4int(header[
"nt"]);
336 header[
"xmin"] * CLHEP::cm, header[
"xmax"] * CLHEP::cm,
337 header[
"ymin"] * CLHEP::cm, header[
"ymax"] * CLHEP::cm,
338 header[
"zmin"] * CLHEP::cm, header[
"zmax"] * CLHEP::cm,
339 header[
"tmin"] * CLHEP::s, header[
"tmax"] * CLHEP::s);
351 if (loopOrder ==
"tzyx")
353 for (G4int i = 0; i < n1; i++)
355 for (G4int j = 0; j < n2; j++)
357 for (G4int k = 0; k < n3; k++)
359 for (G4int l = 0; l < n4; l++)
361 if (!std::getline(file, line))
362 {Terminate(functionName +
"unexpected end to file on line number " + std::to_string(currentLineNumber));}
364 if (std::all_of(line.begin(), line.end(), isspace))
366 ProcessData(line, xIndex, yIndex, zIndex);
367 (*result)(i, j, k, l) = fv;
368 float mag = fv.mag();
369 maximumFieldValue = std::max(maximumFieldValue, mag);
370 minimumFieldValue = std::min(minimumFieldValue, mag);
378 for (G4int l = 0; l < n4; l++)
380 for (G4int k = 0; k < n3; k++)
382 for (G4int j = 0; j < n2; j++)
384 for (G4int i = 0; i < n1; i++)
386 if (!std::getline(file, line))
387 {Terminate(functionName +
"unexpected end to file on line number " + std::to_string(currentLineNumber));}
389 if (std::all_of(line.begin(), line.end(), isspace))
391 ProcessData(line, xIndex, yIndex, zIndex);
392 (*result)(i, j, k, l) = fv;
393 float mag = fv.mag();
394 maximumFieldValue = std::max(maximumFieldValue, mag);
395 minimumFieldValue = std::min(minimumFieldValue, mag);
403 G4cout << functionName <<
"Loaded " << currentLineNumber <<
" lines from file" << G4endl;
404 G4cout << functionName <<
"(Min | Max) field magnitudes in loaded file (before scaling): (" << minimumFieldValue <<
" | " << maximumFieldValue <<
")" << G4endl;
409 const unsigned long xIndex,
410 const unsigned long yIndex,
411 const unsigned long zIndex)
413 std::istringstream liness(line);
415 std::fill(lineData.begin(), lineData.end(), 0);
418 for (
unsigned long i = 1; i < nColumns+1; ++i)
421 if (i < indexOfFirstFieldValue)
422 {value *= CLHEP::cm;}
1D array with spatial mapping derived from BDSArray4DCoords.
2D array with spatial mapping derived from BDSArray4DCoords.
3D array with spatial mapping derived from BDSArray4DCoords.
Overlay of 4D array that provides uniform only spatial coordinate mapping.
General exception with possible name of object and message.
Loader for BDSIM format fields.
BDSArray3DCoords * Load3D(const G4String &fileName)
Load a 3D array.
void CleanUp()
Ensure any member variables are reset between usages.
BDSArray4DCoords * Load4D(const G4String &fileName)
Load a 4D array.
BDSArray1DCoords * Load1D(const G4String &fileName)
Load a 1D array.
void Terminate(const G4String &message="")
Close file and exit program in case of an error.
void Load(const G4String &fileName, const unsigned int nDim)
General loader for any number of dimensions.
void ProcessData(const std::string &line, const unsigned long xIndex, const unsigned long yIndex=0, const unsigned long zIndex=0)
BDSArray2DCoords * Load2D(const G4String &fileName)
Load a 2D array.
BDSDimensionType DetermineDimensionType(G4String dimensionType)
Determine the output format to be used from the input string.
G4String LowerCase(const G4String &str)
Utility function to simplify lots of syntax changes for pedantic g4 changes.