Skip to content

Commit fc23a36

Browse files
maureeungaromaureeungaro
committed
added field inquiry option
1 parent 1da2b31 commit fc23a36

7 files changed

Lines changed: 220 additions & 20 deletions

File tree

‎gemc.cc‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include "gaction.h"
1818
#include "gstreamer.h"
1919
#include "gsystem_options.h"
20+
#include "gfield_options.h"
2021

2122
// c++
2223
#include <cstdio>
@@ -29,11 +30,12 @@ int main(int argc, char* argv[]) {
2930
auto gopts = std::make_shared<GOptions>(argc, argv, base_schema);
3031
auto log = std::make_shared<GLogger>(gopts, SFUNCTION_NAME, GENERAL_LOGGER);
3132

33+
if (gfields::runFieldQueries(gopts)) { return EXIT_SUCCESS; }
34+
3235
auto gui = gopts->getSwitch("gui");
3336
auto nthreads = gemc::get_nthreads(gopts, log);
3437
auto has_startup_geometry = !gsystem::getSystems(gopts).empty();
3538

36-
3739
// createQtApplication returns a QApplication if gui is true
3840
// otherwise, it returns a QCoreApplication and sets the Geant4 CoutDestination to a GBatch_Session
3941
auto app = gemc::makeQtApplication(argc, argv, gui);
@@ -46,7 +48,8 @@ int main(int argc, char* argv[]) {
4648
gemc::start_random_engine(gopts, log);
4749

4850
// init geant4 run manager with then number of threads coming from options. always fails if unavailable
49-
auto runManager = std::unique_ptr<G4RunManager>(G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default, true, nthreads));
51+
auto runManager = std::unique_ptr<G4RunManager>(
52+
G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default, true, nthreads));
5053

5154
// Pre-load streamer plugins before Geant4 creates worker threads. Sanitized Linux
5255
// builds can fail late dlopen() calls from workers with static TLS exhaustion.

‎gemc/gfields/examples/test_gfield_dipole.cc‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,8 @@ int main(int argc, char* argv[]) {
5353
// Initialize GOptions (parsed from YAML or another source)
5454
auto gopts = std::make_shared<GOptions>(argc, argv, gfields::defineOptions());
5555

56+
if (gfields::runFieldQueries(gopts)) { return EXIT_SUCCESS; }
57+
5658
// Create a GMagneto instance to manage fields
5759
auto magneto = std::make_shared<GMagneto>(gopts);
5860

‎gemc/gfields/gfield_options.cc‎

Lines changed: 170 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,116 @@
11
#include "gfield_options.h"
22
#include "gfieldConventions.h"
3+
#include "gmagneto.h"
4+
5+
// geant4
6+
#include "G4UnitsTable.hh"
37

48
// gemc
59
#include "gutilities.h"
610
#include "gfactory_options.h"
11+
#include "goptionsConventions.h"
12+
13+
// c++
14+
#include <cmath>
15+
#include <cstdlib>
16+
#include <fstream>
17+
#include <iostream>
18+
#include <sstream>
719

8-
// namespace to define options
920
namespace gfields {
1021

11-
// Build field definitions by reading the option tree and translating each entry into a GFieldDefinition (non-Doxygen summary).
22+
namespace {
23+
24+
struct FieldQueryPoint {
25+
double position[3] = {0.0, 0.0, 0.0};
26+
std::string source;
27+
int line = 0;
28+
};
29+
30+
bool is_query_set(const std::string& value) {
31+
return !value.empty() && value != UNINITIALIZEDSTRINGQUANTITY && value != "not provided";
32+
}
33+
34+
FieldQueryPoint parse_field_query_point(const std::string& expression, const std::string& source, int line) {
35+
std::string cleaned = expression;
36+
if (const auto comment = cleaned.find('#'); comment != std::string::npos) {
37+
cleaned = cleaned.substr(0, comment);
38+
}
39+
for (auto& c : cleaned) {
40+
if (c == ',') { c = ' '; }
41+
}
42+
cleaned = gutilities::removeLeadingAndTrailingSpacesFromString(cleaned);
43+
44+
FieldQueryPoint point;
45+
point.source = source;
46+
point.line = line;
47+
48+
if (cleaned.empty()) { return point; }
49+
50+
std::istringstream tokens(cleaned);
51+
std::string x;
52+
std::string y;
53+
std::string z;
54+
std::string extra;
55+
tokens >> x >> y >> z >> extra;
56+
if (x.empty() || y.empty() || z.empty() || !extra.empty()) {
57+
std::cerr << FATALERRORL << "field query point must contain exactly three coordinates with units";
58+
if (line > 0) { std::cerr << " at " << source << ":" << line; }
59+
else { std::cerr << " in " << source; }
60+
std::cerr << ". Got <" << expression << ">." << std::endl;
61+
std::exit(EC__G4NUMBERERROR);
62+
}
63+
64+
point.position[0] = gutilities::getG4Number(x, true);
65+
point.position[1] = gutilities::getG4Number(y, true);
66+
point.position[2] = gutilities::getG4Number(z, true);
67+
return point;
68+
}
69+
70+
bool is_blank_or_comment_line(const std::string& line) {
71+
const auto trimmed = gutilities::removeLeadingAndTrailingSpacesFromString(line);
72+
return trimmed.empty() || trimmed[0] == '#';
73+
}
74+
75+
void append_field_query_file_points(const std::string& filename, std::vector<FieldQueryPoint>& points) {
76+
std::ifstream input(filename);
77+
if (!input) {
78+
std::cerr << FATALERRORL << "can't open field query point file " << filename << "." << std::endl;
79+
std::exit(EC__FILENOTFOUND);
80+
}
81+
82+
std::string line;
83+
int line_number = 0;
84+
while (std::getline(input, line)) {
85+
++line_number;
86+
if (is_blank_or_comment_line(line)) { continue; }
87+
points.push_back(parse_field_query_point(line, filename, line_number));
88+
}
89+
}
90+
91+
void print_field_query_result(const std::string& field_name,
92+
const FieldQueryPoint& point,
93+
const double bfield[3]) {
94+
const double bmag = std::sqrt(
95+
bfield[0] * bfield[0] +
96+
bfield[1] * bfield[1] +
97+
bfield[2] * bfield[2]);
98+
std::cout << "field=" << field_name
99+
<< " source=" << point.source;
100+
if (point.line > 0) { std::cout << ":" << point.line; }
101+
std::cout << " x=" << G4BestUnit(point.position[0], "Length")
102+
<< " y=" << G4BestUnit(point.position[1], "Length")
103+
<< " z=" << G4BestUnit(point.position[2], "Length")
104+
<< " Bx=" << G4BestUnit(bfield[0], "Magnetic flux density")
105+
<< " By=" << G4BestUnit(bfield[1], "Magnetic flux density")
106+
<< " Bz=" << G4BestUnit(bfield[2], "Magnetic flux density")
107+
<< " |B|=" << G4BestUnit(bmag, "Magnetic flux density")
108+
<< std::endl;
109+
}
110+
111+
} // namespace
112+
113+
// Build field definitions by reading the option tree and translating each entry into a GFieldDefinition.
12114
std::vector<GFieldDefinition> get_GFieldDefinition(const std::shared_ptr<GOptions>& gopts) {
13115
std::vector<GFieldDefinition> gfield_defs;
14116

@@ -19,20 +121,31 @@ std::vector<GFieldDefinition> get_GFieldDefinition(const std::shared_ptr<GOption
19121
GFieldDefinition gfield_def = GFieldDefinition();
20122

21123
// Core identity and integration configuration.
22-
gfield_def.name = gopts->get_variable_in_option<std::string>(gmultipoles_item, "name", goptions::NODFLT);
23-
gfield_def.integration_stepper = gopts->get_variable_in_option<std::string>(gmultipoles_item, "integration_stepper", GFIELD_DEFAULT_INTEGRATION_STEPPER);
24-
gfield_def.minimum_step = gutilities::getG4Number(gopts->get_variable_in_option<std::string>(gmultipoles_item, "minimum_step", GFIELD_DEFAULT_MINIMUM_STEP));
124+
gfield_def.name = gopts->get_variable_in_option<std::string>(
125+
gmultipoles_item, "name", goptions::NODFLT);
126+
gfield_def.integration_stepper = gopts->get_variable_in_option<std::string>(
127+
gmultipoles_item, "integration_stepper", GFIELD_DEFAULT_INTEGRATION_STEPPER);
128+
gfield_def.minimum_step = gutilities::getG4Number(gopts->get_variable_in_option<std::string>(
129+
gmultipoles_item, "minimum_step", GFIELD_DEFAULT_MINIMUM_STEP));
25130

26131
// Multipole parameters:
27132
// Values are stored as strings to preserve unit expressions and are parsed later by the concrete field.
28-
gfield_def.add_map_parameter("pole_number", gopts->get_variable_in_option<std::string>(gmultipoles_item, "pole_number", goptions::NODFLT));
29-
gfield_def.add_map_parameter("vx", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vx", GFIELD_DEFAULT_VERTEX));
30-
gfield_def.add_map_parameter("vy", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vy", GFIELD_DEFAULT_VERTEX));
31-
gfield_def.add_map_parameter("vz", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vz", GFIELD_DEFAULT_VERTEX));
32-
gfield_def.add_map_parameter("rotation_angle", gopts->get_variable_in_option<std::string>(gmultipoles_item, "rotation_angle", GFIELD_DEFAULT_ROTANGLE));
33-
gfield_def.add_map_parameter("rotaxis", gopts->get_variable_in_option<std::string>(gmultipoles_item, "rotaxis", goptions::NODFLT));
34-
gfield_def.add_map_parameter("strength", gopts->get_variable_in_option<std::string>(gmultipoles_item, "strength", goptions::NODFLT));
35-
gfield_def.add_map_parameter("longitudinal", gopts->get_variable_in_option<std::string>(gmultipoles_item, "longitudinal", "false"));
133+
gfield_def.add_map_parameter("pole_number", gopts->get_variable_in_option<std::string>(
134+
gmultipoles_item, "pole_number", goptions::NODFLT));
135+
gfield_def.add_map_parameter("vx", gopts->get_variable_in_option<std::string>(
136+
gmultipoles_item, "vx", GFIELD_DEFAULT_VERTEX));
137+
gfield_def.add_map_parameter("vy", gopts->get_variable_in_option<std::string>(
138+
gmultipoles_item, "vy", GFIELD_DEFAULT_VERTEX));
139+
gfield_def.add_map_parameter("vz", gopts->get_variable_in_option<std::string>(
140+
gmultipoles_item, "vz", GFIELD_DEFAULT_VERTEX));
141+
gfield_def.add_map_parameter("rotation_angle", gopts->get_variable_in_option<std::string>(
142+
gmultipoles_item, "rotation_angle", GFIELD_DEFAULT_ROTANGLE));
143+
gfield_def.add_map_parameter("rotaxis", gopts->get_variable_in_option<std::string>(
144+
gmultipoles_item, "rotaxis", goptions::NODFLT));
145+
gfield_def.add_map_parameter("strength", gopts->get_variable_in_option<std::string>(
146+
gmultipoles_item, "strength", goptions::NODFLT));
147+
gfield_def.add_map_parameter("longitudinal", gopts->get_variable_in_option<std::string>(
148+
gmultipoles_item, "longitudinal", "false"));
36149

37150
// The type field controls the shared-library plugin name through GFieldDefinition::gfieldPluginName().
38151
gfield_def.type = "multipoles";
@@ -71,6 +184,50 @@ GOptions defineOptions() {
71184
};
72185
goptions.defineOption("gmultipoles", "define the e.m. gmultipoles", gmultipoles, help);
73186

187+
goptions.defineOption(
188+
GVariable("fieldAt", UNINITIALIZEDSTRINGQUANTITY, "query all configured fields at x y z"),
189+
"Evaluate all configured electromagnetic fields at one absolute coordinate.\n \n"
190+
"The value must contain three coordinate expressions with units, separated by spaces.\n \n"
191+
"Example: -fieldAt=\"10*cm 0*mm 2*m\"\n \n");
192+
193+
goptions.defineOption(
194+
GVariable("fieldMapPoints", UNINITIALIZEDSTRINGQUANTITY, "ASCII file of x y z points for field queries"),
195+
"Evaluate all configured electromagnetic fields at coordinates listed in an ASCII file.\n \n"
196+
"Each non-empty, non-comment line must contain three coordinate expressions with units.\n"
197+
"Coordinates may be separated by spaces or commas. Lines beginning with # are ignored.\n \n"
198+
"Example: -fieldMapPoints=points.txt\n \n");
199+
74200
return goptions;
75201
}
202+
203+
bool runFieldQueries(const std::shared_ptr<GOptions>& gopts) {
204+
const auto field_at = gopts->getScalarString("fieldAt");
205+
const auto field_map_points = gopts->getScalarString("fieldMapPoints");
206+
207+
if (!is_query_set(field_at) && !is_query_set(field_map_points)) { return false; }
208+
209+
std::vector<FieldQueryPoint> points;
210+
if (is_query_set(field_at)) { points.push_back(parse_field_query_point(field_at, "fieldAt", 0)); }
211+
if (is_query_set(field_map_points)) { append_field_query_file_points(field_map_points, points); }
212+
213+
auto magneto = std::make_shared<GMagneto>(gopts);
214+
auto field_names = magneto->getFieldNames();
215+
if (field_names.empty()) {
216+
std::cerr << FATALERRORL << "field query requested, but no electromagnetic fields are configured."
217+
<< std::endl;
218+
std::exit(EC__NOOPTIONFOUND);
219+
}
220+
221+
std::cout << "# field query results" << std::endl;
222+
for (const auto& point : points) {
223+
for (const auto& field_name : field_names) {
224+
double bfield[3] = {0.0, 0.0, 0.0};
225+
magneto->getField(field_name)->GetFieldValue(point.position, bfield);
226+
print_field_query_result(field_name, point, bfield);
227+
}
228+
}
229+
230+
return true;
231+
}
232+
76233
}

‎gemc/gfields/gfield_options.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,17 @@ std::vector<GFieldDefinition> get_GFieldDefinition(const std::shared_ptr<GOption
3232
*/
3333
GOptions defineOptions();
3434

35+
/**
36+
* \brief Evaluate configured fields from \c fieldAt and \c fieldMapPoints options, if requested.
37+
* \param gopts Shared options container with field definitions and query options.
38+
* \return True when at least one query option was present and processed.
39+
*
40+
* The query options are intended for command-line inspection workflows:
41+
* - \c fieldAt accepts one coordinate triplet with units, for example \c "10*cm 0*mm 2*m".
42+
* - \c fieldMapPoints accepts an ASCII file with one \c x y z coordinate triplet per line.
43+
*/
44+
bool runFieldQueries(const std::shared_ptr<GOptions>& gopts);
45+
3546
} // namespace gfields
3647

3748
/** @} */

‎gemc/gfields/gmagneto.cc‎

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,9 @@
55
#include "gmagneto.h"
66
#include "gfield_options.h"
77

8+
// c++
9+
#include <algorithm>
10+
811
// #include "G4TransportationManager.hh"
912
// #include "G4PropagatorInField.hh"
1013

@@ -28,7 +31,8 @@ GMagneto::GMagneto(const std::shared_ptr<GOptions>& gopts) : GBase(gopts, GMAGNE
2831
// Only create each named field once; repeated names are ignored by this map check.
2932
if (fields_map->find(name) == fields_map->end()) {
3033
// Load the plugin, instantiate the field object, and cache it by name.
31-
fields_map->emplace(name, gFieldManager.LoadAndRegisterObjectFromLibrary<GField>(field_definition.gfieldPluginName(), gopts));
34+
fields_map->emplace(name, gFieldManager.LoadAndRegisterObjectFromLibrary<GField>(
35+
field_definition.gfieldPluginName(), gopts));
3236

3337
// Pass the configuration down to the concrete implementation so it can parse/cache parameters.
3438
fields_map->at(name)->load_field_definitions(field_definition);
@@ -41,3 +45,11 @@ GMagneto::GMagneto(const std::shared_ptr<GOptions>& gopts) : GBase(gopts, GMAGNE
4145
// TODO: add min and max steps
4246
// G4TransportationManager::GetTransportationManager()->GetPropagatorInField()->SetLargestAcceptableStep(10);
4347
}
48+
49+
std::vector<std::string> GMagneto::getFieldNames() const {
50+
std::vector<std::string> names;
51+
names.reserve(fields_map->size());
52+
for (const auto& [name, field] : *fields_map) { names.push_back(name); }
53+
std::sort(names.begin(), names.end());
54+
return names;
55+
}

‎gemc/gfields/gmagneto.h‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -26,11 +26,12 @@
2626
*
2727
* Ownership model:
2828
* - \ref GMagneto "GMagneto" owns a map of field objects (\ref GField) created via dynamic plugin loading.
29-
* - For each field it also owns a corresponding \c G4FieldManager created by \ref GField::create_FieldManager "create_FieldManager()".
29+
* - For each field it also owns a corresponding \c G4FieldManager created by
30+
* \ref GField::create_FieldManager "create_FieldManager()".
3031
*
3132
* Lifecycle:
32-
* - Fields and managers are constructed during \ref GMagneto::GMagneto "GMagneto()" based on the field definitions
33-
* produced by gfields::get_GFieldDefinition().
33+
* - Fields and managers are constructed during \ref GMagneto::GMagneto "GMagneto()" based on the field
34+
* definitions produced by gfields::get_GFieldDefinition().
3435
* - Maps live for the lifetime of the \ref GMagneto "GMagneto" instance.
3536
*/
3637
class GMagneto : public GBase<GMagneto> {
@@ -65,6 +66,12 @@ class GMagneto : public GBase<GMagneto> {
6566
*/
6667
bool isField(const std::string& name) const { return fields_map->find(name) != fields_map->end(); }
6768

69+
/**
70+
* \brief Return the configured field names.
71+
* \return Sorted list of field names loaded by this manager.
72+
*/
73+
std::vector<std::string> getFieldNames() const;
74+
6875
/**
6976
* \brief Retrieve a field object by name.
7077
* \param name Field name key.
@@ -93,7 +100,9 @@ class GMagneto : public GBase<GMagneto> {
93100
log->info(2, "GFieldManager >", key, "< >", value);
94101
}
95102

96-
if (fields_manager->find(name) == fields_manager->end()) { log->error(ERR_WRONG_FIELD_NOT_FOUND, "GField >", name, "< not found. Exiting."); }
103+
if (fields_manager->find(name) == fields_manager->end()) {
104+
log->error(ERR_WRONG_FIELD_NOT_FOUND, "GField >", name, "< not found. Exiting.");
105+
}
97106

98107
return fields_manager->at(name);
99108
}

‎gemc/gfields/meson.build‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ verbosities = ['-verbosity.gfield=2',
66
'-debug.gfield=true'
77
]
88
input_file = [meson.current_source_dir() + '/examples/dipole.yaml']
9+
field_points_file = meson.current_source_dir() + '/examples/field_points.txt'
910

1011

1112
multipoles_files = [files(
@@ -38,5 +39,10 @@ LD += {
3839
'examples' : {
3940
'test_gfield_dipole' : [example_source, input_file],
4041
'test_gfield_dipole_verbose' : [example_source, input_file + verbosities],
42+
'test_gfield_dipole_field_at' : [example_source, input_file + ['-fieldAt=10*cm 0*cm 0*cm']],
43+
'test_gfield_dipole_field_map_points' : [
44+
example_source,
45+
input_file + ['-fieldMapPoints=' + field_points_file],
46+
],
4147
}
4248
}

0 commit comments

Comments
 (0)