Skip to content

Commit 19bbfbe

Browse files
maureeungaromaureeungaro
committed
fixed gfield multipoles and added longitudinal option. added ascii_db option to define additional search path for ascii. using assimp internal zlib for non macos systems
1 parent 1ca3234 commit 19bbfbe

18 files changed

Lines changed: 223 additions & 130 deletions

File tree

‎ci/build.sh‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,10 +5,16 @@
55
# Container run:
66
# docker run -it --rm --platform linux/amd64 jeffersonlab/geant4:g4v11.3.2-almalinux94 sh
77
# git clone http://github.com/gemc/src /root/src && cd /root/src
8-
# ./ci/build.sh none
8+
# ./ci/build.sh
99

1010
source ci/env.sh
1111

12+
# module gemc gives $GEMC (used in meson prefix) and PKG_CONFIG_PATH
13+
# not strickly necessary, one could set those manually
14+
module load gemc/dev3
15+
echo GEMC for prefix set to: $GEMC
16+
echo
17+
1218
meson_option=$(meson_options $1)
1319
max_threads=$(max_j)
1420

‎ci/env.sh‎

Lines changed: 0 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -87,19 +87,6 @@ else
8787
echo "GITHUB_WORKFLOW: ${GITHUB_WORKFLOW}"
8888
fi
8989
source /etc/profile.d/localSetup.sh
90-
module load gemc/dev3
9190
echo
9291
fi
9392

94-
mkdir -p $GEMC/lib/pkgconfig
95-
96-
# if $GEMC/lib/pkgconfig/geant4.pc is not existing, run meson/install_geant4_root_pkgconfig.py
97-
if [ ! -f $GEMC/lib/pkgconfig/geant4.pc ]; then
98-
echo " > Running meson/install_geant4_root_pkgconfig.py"
99-
python3 meson/install_geant4_root_pkgconfig.py
100-
fi
101-
102-
# since the pkgconfig files are installed after the module loads, we need to reload the modules
103-
module unload gemc
104-
module load gemc/dev3
105-

‎examples/geant4_basic/b2/b2.yaml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@ gmultipoles:
2828
- name: dipole
2929
pole_number: 2
3030
strength: 2
31+
rotaxis: z
3132
rotation_angle: "30*deg"
3233

3334
root: G4Box, 12*cm, 12*cm, 18*cm, G4_WATER

‎gfields/examples/dipole.yaml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,5 @@ gmultipoles:
22
- name: dipole
33
pole_number: 2
44
strength: 2
5+
rotaxis: z
56
rotation_angle: "30*deg"

‎gfields/gfieldConventions.h‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,11 +7,8 @@
77
// other default parameters are in the form of strings
88

99
// multipoles.
10-
#define GFIELD_DEFAULT_POLE_NUMBER "1"
1110
#define GFIELD_DEFAULT_VERTEX "0*mm"
1211
#define GFIELD_DEFAULT_ROTANGLE "0*deg"
13-
#define GFIELD_DEFAULT_ROTAXIS "Z"
14-
#define GFIELD_DEFAULT_STRENGTH "1.0"
1512
#define GFIELD_DEFAULT_VERBOSITY "0"
1613

1714
// error codes in the 1200

‎gfields/gfieldFactories/multipoles/gfield_multipoles.cc‎

Lines changed: 133 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -9,88 +9,154 @@
99
#include "gfieldConventions.h"
1010

1111
// gemv
12-
#include "gutilities.h"
12+
// #include "gutilities.h"
1313

14-
// c++
15-
#include <iostream>
16-
17-
using namespace std;
1814

1915
// tells the DLL how to create a GFieldFactory
20-
extern "C" GField *GFieldFactory(void) {
21-
return static_cast<GField *>(new GField_MultipolesFactory);
22-
}
16+
extern "C" GField* GFieldFactory(void) { return static_cast<GField*>(new GField_MultipolesFactory); }
2317

2418

2519
// for now this implementation follows gemc
26-
// reference of this implementation: https://uspas.fnal.gov/materials/12MSU/magnet_elements.pdf
27-
void GField_MultipolesFactory::GetFieldValue(const double pos[3], G4double *bfield) const {
28-
29-
G4ThreeVector x0(pos[0], pos[1], pos[2]);
30-
G4ThreeVector x1(origin[0], origin[1], origin[2]);
31-
G4ThreeVector x2;
32-
G4ThreeVector x0_local;
33-
if (rotaxis == 0) {
34-
x2 = G4ThreeVector(0 * CLHEP::cm, 0 * CLHEP::cm, 1 * CLHEP::cm).rotateX(rotation_angle) + x1;
35-
x0_local = (x0 - x1).rotateX(-rotation_angle);
36-
} else if (rotaxis == 1) {
37-
x2 = G4ThreeVector(0 * CLHEP::cm, 0 * CLHEP::cm, 1 * CLHEP::cm).rotateY(rotation_angle) + x1;
38-
x0_local = (x0 - x1).rotateY(-rotation_angle);
39-
} else if (rotaxis == 2) {
40-
x2 = G4ThreeVector(0 * CLHEP::cm, 0 * CLHEP::cm, 1 * CLHEP::cm).rotateZ(rotation_angle) + x1;
41-
x0_local = (x0 - x1).rotateZ(-rotation_angle);
20+
// references of this implementation:
21+
// - https://cds.cern.ch/record/1333874/files/1.pdf
22+
// - https://uspas.fnal.gov/materials/12MSU/magnet_elements.pdf
23+
// - https://cas.web.cern.ch/sites/default/files/lectures/bruges-2009/wolski-1.pdf
24+
// notice strength is defined at a reference radius of 1m
25+
void GField_MultipolesFactory::GetFieldValue(const double pos[3], G4double* bfield) const
26+
{
27+
// ======= Configuration / conventions =======
28+
// strength: Tesla at reference radius r0 for all multipole orders
29+
const double r0 = CLHEP::m; // reference radius; make this a member if you want
30+
31+
// ======= Basic checks =======
32+
if (pole_number < 2 || (pole_number % 2) != 0) {
33+
log->error(ERR_WRONG_POLE_NUMBER,
34+
"Pole number must be an even integer >= 2 (2=dipole,4=quadrupole,...)");
35+
bfield[0] = bfield[1] = bfield[2] = 0.0;
36+
return;
37+
}
38+
39+
// ======= Positions and local frame =======
40+
const G4ThreeVector x0(pos[0], pos[1], pos[2]); // query point (lab)
41+
const G4ThreeVector x1(origin[0], origin[1], origin[2]); // magnet origin (lab)
42+
43+
// shift to magnet-centered coordinates and "unroll" the magnet by -rotation_angle
44+
G4ThreeVector p = x0 - x1;
45+
if (rotaxis == 0) p.rotateX(-rotation_angle);
46+
else if (rotaxis == 1) p.rotateY(-rotation_angle);
47+
else if (rotaxis == 2) p.rotateZ(-rotation_angle);
48+
49+
// Identify transverse plane (u,v) ⟂ to the axis, and the axial coordinate w (unused here)
50+
double u=0.0, v=0.0;
51+
switch (rotaxis) {
52+
case 0: u = p.y(); v = p.z(); break; // axis = X ⇒ transverse plane = (Y,Z)
53+
case 1: u = p.z(); v = p.x(); break; // axis = Y ⇒ transverse plane = (Z,X)
54+
case 2: u = p.x(); v = p.y(); break; // axis = Z ⇒ transverse plane = (X,Y)
55+
}
56+
57+
const double r = std::hypot(u, v); // transverse radius
58+
const double phi = std::atan2(v, u); // azimuth in transverse plane
59+
60+
// ======= Axial (solenoid-like) mode if explicitly requested =======
61+
if (longitudinal) {
62+
// Uniform axial field aligned with rotaxis; not a multipole.
63+
G4ThreeVector B_local(0,0,0);
64+
switch (rotaxis) {
65+
case 0: B_local.setX(strength); break;
66+
case 1: B_local.setY(strength); break;
67+
case 2: B_local.setZ(strength); break;
68+
}
69+
// Roll back to lab
70+
G4ThreeVector B_lab = B_local;
71+
if (rotaxis == 0) B_lab.rotateX(+rotation_angle);
72+
else if (rotaxis == 1) B_lab.rotateY(+rotation_angle);
73+
else if (rotaxis == 2) B_lab.rotateZ(+rotation_angle);
74+
75+
bfield[0] = B_lab.x();
76+
bfield[1] = B_lab.y();
77+
bfield[2] = B_lab.z();
78+
79+
log->info(2, "Axial field mode (solenoid-like). Strength: ", strength,
80+
" T, Field: (", bfield[0], ", ", bfield[1], ", ", bfield[2], ")");
81+
return;
82+
}
83+
84+
// ======= Transverse multipole (standard accelerator definition) =======
85+
const int n = pole_number / 2; // 1=dipole, 2=quadrupole, 3=sextupole, ...
86+
const int a = n - 1; // power of r
87+
88+
// r^a scaling with reference radius r0; for dipole (a=0) this is 1
89+
double ra = 1.0;
90+
if (a > 0) {
91+
// If r=0 and a>0, |B|=0 (for ideal multipoles).
92+
if (r == 0.0) {
93+
bfield[0] = bfield[1] = bfield[2] = 0.0;
94+
return;
95+
}
96+
ra = std::pow(r / r0, static_cast<double>(a));
4297
}
4398

44-
G4double r = (x2 - x1).cross(x1 - x0).mag() / (x2 - x1).mag(); //distance from x0 to line x1-x2
45-
G4double phi = atan2(x0_local.y(), x0_local.x());
46-
47-
G4ThreeVector B_local;
48-
if (pole_number == 2) {
49-
B_local.setX(0);
50-
B_local.setY(strength);
51-
B_local.setZ(0);
52-
} else if (pole_number > 0) {
53-
int a = pole_number / 2 - 1;
54-
B_local.setX(strength * pow(r / CLHEP::m, a) * sin(a * phi));
55-
B_local.setY(strength * pow(r / CLHEP::m, a) * cos(a * phi));
56-
B_local.setZ(0);
57-
} else {
58-
log->error(ERR_WRONG_POLE_NUMBER, "GField_MultipolesFactory::GetFieldValue: Pole number " + to_string(pole_number) + " not supported. Exiting.");
99+
// "Normal" multipole angular dependence:
100+
// Bu = strength * r^(n-1) * cos((n-1)*phi)
101+
// Bv = strength * r^(n-1) * sin((n-1)*phi)
102+
// This yields a constant transverse dipole for n=1.
103+
const double Bu = strength * ra * std::cos(a * phi);
104+
const double Bv = strength * ra * std::sin(a * phi);
105+
106+
// Place (Bu,Bv) into the correct transverse components of B_local; axial component = 0
107+
G4ThreeVector B_local(0,0,0);
108+
switch (rotaxis) {
109+
case 0: B_local.setY(Bu); B_local.setZ(Bv); break; // axis X → (Y,Z)
110+
case 1: B_local.setZ(Bu); B_local.setX(Bv); break; // axis Y → (Z,X)
111+
case 2: B_local.setX(Bu); B_local.setY(Bv); break; // axis Z → (X,Y)
59112
}
60113

114+
// Rotate (roll) back to lab
61115
G4ThreeVector B_lab = B_local;
62-
if (rotaxis == 0) { B_lab.rotateX(rotation_angle); }
63-
else if (rotaxis == 1) { B_lab.rotateY(rotation_angle); }
64-
else if (rotaxis == 2) { B_lab.rotateZ(rotation_angle); }
116+
if (rotaxis == 0) B_lab.rotateX(+rotation_angle);
117+
else if (rotaxis == 1) B_lab.rotateY(+rotation_angle);
118+
else if (rotaxis == 2) B_lab.rotateZ(+rotation_angle);
65119

66-
bfield[0] = B_lab.x();
67-
bfield[1] = B_lab.y();
68-
bfield[2] = B_lab.z();
120+
// Output
121+
bfield[0] = B_lab.x();
122+
bfield[1] = B_lab.y();
123+
bfield[2] = B_lab.z();
69124

70125
log->info(2, "Pole Number: ", pole_number,
71-
", Strength: ", strength,
72-
", Requested at: (", pos[0], ", ", pos[1], ", ", pos[2], ")",
73-
", Rotation angle: ", rotation_angle,
74-
", Rotation axis: ", rotaxis,
75-
", Field: (", bfield[0], ", ", bfield[1], ", ", bfield[2], ")");
126+
", n: ", n,
127+
", Strength: ", strength,
128+
", Requested at: (", pos[0], ", ", pos[1], ", ", pos[2], ")",
129+
", Rotation angle: ", rotation_angle,
130+
", Rotation axis: ", rotaxis,
131+
", longitudinal: ", longitudinal,
132+
", Field: (", bfield[0], ", ", bfield[1], ", ", bfield[2], ")");
76133
}
77134

135+
78136
void GField_MultipolesFactory::load_field_definitions(GFieldDefinition gfd) {
79-
gfield_definitions = gfd;
80-
81-
pole_number = get_field_parameter_int("pole_number");
82-
origin[0] = get_field_parameter_double("vx");
83-
origin[1] = get_field_parameter_double("vy");
84-
origin[2] = get_field_parameter_double("vz");
85-
rotation_angle = get_field_parameter_double("rotation_angle");
86-
if( gfield_definitions.field_parameters["rotaxis"] == "X") {
87-
rotaxis = 0;
88-
} else if( gfield_definitions.field_parameters["rotaxis"] == "Y") {
89-
rotaxis = 1;
90-
} else if( gfield_definitions.field_parameters["rotaxis"] == "Z") {
91-
rotaxis = 2;
92-
} else {
93-
log->error(ERR_WRONG_FIELD_ROTATION, "GField_MultipolesFactory::load_field_definitions: Rotation axis " + gfield_definitions.field_parameters["rotaxis"] + " not supported. Exiting.");
94-
}
95-
strength = get_field_parameter_double("strength");
137+
gfield_definitions = gfd;
138+
139+
pole_number = get_field_parameter_int("pole_number");
140+
origin[0] = get_field_parameter_double("vx");
141+
origin[1] = get_field_parameter_double("vy");
142+
origin[2] = get_field_parameter_double("vz");
143+
rotation_angle = get_field_parameter_double("rotation_angle");
144+
auto rot_axis_option = gfield_definitions.field_parameters["rotaxis"];
145+
longitudinal = false;
146+
if ( gfield_definitions.field_parameters["longitudinal"] == "true" ) {
147+
longitudinal = true;
148+
log->info(1, "Longitudinal field");
149+
} else {
150+
log->info(1, "Transverse field");
151+
}
152+
153+
if (rot_axis_option == "X" || rot_axis_option == "x") { rotaxis = 0; }
154+
else if (rot_axis_option == "Y" || rot_axis_option == "y") { rotaxis = 1; }
155+
else if (rot_axis_option == "Z" || rot_axis_option == "z") { rotaxis = 2; }
156+
else {
157+
log->error(ERR_WRONG_FIELD_ROTATION, "GField_MultipolesFactory::load_field_definitions: Rotation axis " + gfield_definitions.field_parameters["rotaxis"] +
158+
" not supported. Exiting.");
159+
}
160+
log->info(1, "Rotation axis: ", rotaxis);
161+
strength = get_field_parameter_double("strength");
96162
}

‎gfields/gfieldFactories/multipoles/gfield_multipoles.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@ class GField_MultipolesFactory : public GField {
3333
G4double rotation_angle; ///< Rotation angle of the field in degrees.
3434
int rotaxis; ///< Axis of rotation: 0 for X, 1 for Y, 2 for Z.
3535
G4double strength; ///< Strength of the multipole field.
36+
bool longitudinal; ///< Longitudinal or transverse field.
3637

3738
};
3839

‎gfields/gfield_options.cc‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -20,13 +20,14 @@ std::vector<GFieldDefinition> get_GFieldDefinition(const std::shared_ptr<GOption
2020
gfield_def.minimum_step = gutilities::getG4Number(gopts->get_variable_in_option<std::string>(gmultipoles_item, "minimum_step", GFIELD_DEFAULT_MINIMUM_STEP));
2121

2222
// if present, add remaining multipoles parameters
23-
gfield_def.add_map_parameter("pole_number", gopts->get_variable_in_option<std::string>(gmultipoles_item, "pole_number", GFIELD_DEFAULT_POLE_NUMBER));
23+
gfield_def.add_map_parameter("pole_number", gopts->get_variable_in_option<std::string>(gmultipoles_item, "pole_number", goptions::NODFLT));
2424
gfield_def.add_map_parameter("vx", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vx", GFIELD_DEFAULT_VERTEX));
2525
gfield_def.add_map_parameter("vy", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vy", GFIELD_DEFAULT_VERTEX));
2626
gfield_def.add_map_parameter("vz", gopts->get_variable_in_option<std::string>(gmultipoles_item, "vz", GFIELD_DEFAULT_VERTEX));
2727
gfield_def.add_map_parameter("rotation_angle", gopts->get_variable_in_option<std::string>(gmultipoles_item, "rotation_angle", GFIELD_DEFAULT_ROTANGLE));
28-
gfield_def.add_map_parameter("rotaxis", gopts->get_variable_in_option<std::string>(gmultipoles_item, "rotaxis", GFIELD_DEFAULT_ROTAXIS));
29-
gfield_def.add_map_parameter("strength", gopts->get_variable_in_option<std::string>(gmultipoles_item, "strength", GFIELD_DEFAULT_STRENGTH));
28+
gfield_def.add_map_parameter("rotaxis", gopts->get_variable_in_option<std::string>(gmultipoles_item, "rotaxis", goptions::NODFLT));
29+
gfield_def.add_map_parameter("strength", gopts->get_variable_in_option<std::string>(gmultipoles_item, "strength", goptions::NODFLT));
30+
gfield_def.add_map_parameter("longitudinal", gopts->get_variable_in_option<std::string>(gmultipoles_item, "longitudinal", "false"));
3031
gfield_def.type = "multipoles";
3132
gfield_defs.push_back(gfield_def);
3233
}
@@ -45,13 +46,14 @@ GOptions defineOptions() {
4546
{"name", goptions::NODFLT, "Field name"},
4647
{"integration_stepper", GFIELD_DEFAULT_INTEGRATION_STEPPER, "Integration stepper"},
4748
{"minimum_step", GFIELD_DEFAULT_MINIMUM_STEP, "Minimum step for the G4ChordFinder"},
48-
{"pole_number", GFIELD_DEFAULT_POLE_NUMBER, "Pole numner"},
49+
{"pole_number", goptions::NODFLT, "Pole numner"},
4950
{"vx", GFIELD_DEFAULT_VERTEX, "x component of the origin of the multipole"},
5051
{"vy", GFIELD_DEFAULT_VERTEX, "y component of the origin of the multipole"},
5152
{"vz", GFIELD_DEFAULT_VERTEX, "z component of the origin of the multipole"},
5253
{"rotation_angle", GFIELD_DEFAULT_ROTANGLE, "rotation angle of the multipole"},
53-
{"rotaxis", GFIELD_DEFAULT_ROTAXIS, "rotation axis of the multipole"},
54-
{"strength", GFIELD_DEFAULT_STRENGTH, "strength of the multipole"}
54+
{"rotaxis", goptions::NODFLT, "rotation axis (roll) of the multipole"},
55+
{"strength", goptions::NODFLT, "strength of the multipole in tesla at radius of 1m"},
56+
{"longitudinal", "false", "if set to true, the field is aligned with rotaxis (solenoid like)"}
5557
};
5658
goptions.defineOption("gmultipoles", "define the e.m. gmultipoles", gmultipoles, help);
5759

‎gfields/gmagneto.cc‎

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

8+
// #include "G4TransportationManager.hh"
9+
// #include "G4PropagatorInField.hh"
10+
811

912
GMagneto::GMagneto(std::shared_ptr<GOptions> gopts) : log(std::make_shared<GLogger>(gopts, GFIELD_LOGGER, "GMagneto")){
1013

@@ -29,4 +32,9 @@ GMagneto::GMagneto(std::shared_ptr<GOptions> gopts) : log(std::make_shared<GLogg
2932
}
3033
}
3134

35+
// TODO: add min and max steps
36+
37+
// G4TransportationManager::GetTransportationManager()->GetPropagatorInField()->SetLargestAcceptableStep(10);
38+
39+
3240
}

‎goptions/goptions.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -586,7 +586,7 @@ void GOptions::print_version() {
586586
cout << endl << asterisks << endl;
587587
cout << " " << KGRN << KBOLD << executableName << RST << " version: " << KGRN << gversion << RST << endl;
588588
cout << " Called from: " << KGRN << executableCallingDir << RST << endl;
589-
cout << " Executed from: " << KGRN << installDir << "/bin" << RST << endl;
589+
cout << " Install: " << KGRN << installDir << "/bin" << RST << endl; //
590590
cout << " Released on: " << KGRN << grelease_date << RST << endl;
591591
cout << " GEMC Reference: " << KGRN << greference << RST << endl;
592592
cout << " GEMC Homepage: " << KGRN << gweb << RST << endl;

0 commit comments

Comments
 (0)