BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSVisFieldModel.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) BDSIM Collaboration, 2001 - 2026.
3
4This file is part of BDSIM.
5
6BDSIM is free software: you can redistribute it and/or modify
7it under the terms of the GNU General Public License as published
8by the Free Software Foundation version 3 of the License.
9
10BDSIM is distributed in the hope that it will be useful, but
11WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
17*/
18#include "BDSColourScaleMagma.hh"
19#include "BDSColourScaleViridis.hh"
20#include "BDSFieldQueryForVis.hh"
21#include "BDSFieldQueryInfo.hh"
22#include "BDSVisFieldModel.hh"
23
24#include "G4ArrowModel.hh"
25#include "G4Box.hh"
26#include "G4Circle.hh"
27#include "G4Colour.hh"
28#include "G4String.hh"
29#include "G4ThreeVector.hh"
30#include "G4Types.hh"
31#include "G4VGraphicsScene.hh"
32#include "G4VMarker.hh"
33#include "G4VisAttributes.hh"
34#include "BDSWarning.hh"
35
36#include <algorithm>
37#include <array>
38#include <cmath>
39#include <cstdio>
40#include <iomanip>
41#include <ios>
42#include <limits>
43#include <sstream>
44#include <string>
45#include <utility>
46#include <vector>
47
48G4int BDSVisFieldModel::instanceCounter = 0;
49
50BDSVisFieldModel::BDSVisFieldModel(const std::vector<BDSFieldQueryInfo*>& queriesIn):
51 queries(queriesIn),
52 pointVisB(nullptr),
53 pointVisE(nullptr)
54{
55 fGlobalTag = "BDSVisFieldModel";
56 fGlobalDescription = "Field view #" + std::to_string(instanceCounter);
57 instanceCounter++;
58
59 auto bFieldColour = BDSColourScaleViridis();
60 pointVisB = new G4VisAttributes(bFieldColour.GetValue(0));
61 auto eFieldColour = BDSColourScaleMagma();
62 pointVisE = new G4VisAttributes(eFieldColour.GetValue(0));
63
64 pointVisB->SetForceSolid(true); // no default argument for g4 V < 10.3, so specify true always
65 pointVisE->SetForceSolid(true);
66}
67
68BDSVisFieldModel::~BDSVisFieldModel()
69{
70 delete pointVisB;
71 delete pointVisE;
72}
73
74void BDSVisFieldModel::DescribeYourselfTo(G4VGraphicsScene& sceneHandler)
75{
76 auto bFieldColour = BDSColourScaleViridis();
77 auto eFieldColour = BDSColourScaleMagma();
78
79 BDSFieldQueryForVis querier;
80 for (const auto& query : queries)
81 {
82 querier.QueryField(query);
83 G4double maxFieldB = querier.MaxFieldB();
84 G4double maxFieldE = querier.MaxFieldE();
85 const std::vector<std::array<G4double, 9>> xyzBEs = querier.Values();
86
87 G4double arrowLength = CalculateArrowLength(query);
88 if (arrowLength > 0.01 * std::numeric_limits<G4double>::max()) // This shouldn't happen, but is a sanity check
89 {arrowLength = 1e-3 * sceneHandler.GetExtent().GetExtentRadius();}
90
91 G4double arrowWidth = 0.3*arrowLength;
92 G4double pointSize = std::max(0.1*arrowWidth, 5.0);
93
94 G4bool drawArrows = query->drawArrows;
95 G4bool drawBoxes = query->drawBoxes;
96
97 if (!drawArrows && !drawBoxes)
98 {
99 BDS::Warning("Neither drawArrows or drawBoxes are turned on in query \"" + query->name + "\" -> skipping");
100 continue; // skip this query
101 }
102
103 G4double boxAlpha = query->boxAlpha;
104 if (boxAlpha <= 0)
105 {
106 BDS::Warning("boxAlpha in query \"" + query->name + "\" must be > 0. Setting to 0.01");
107 boxAlpha = 0.01;
108 }
109 if (boxAlpha > 1)
110 {
111 BDS::Warning("boxAlpha in query \"" + query->name + "\" must be <= 1. Setting to 1.0");
112 boxAlpha = 1.0;
113 }
114
115 G4ThreeVector boxHalfSize = BoxHalfSize(query);
116 G4AffineTransform qTransform = query->globalTransform;
117 G4RotationMatrix qRotation = qTransform.NetRotation().inverse();
118 G4Box aBox("b", boxHalfSize.x(), boxHalfSize.y(), boxHalfSize.z());
119 G4Polyhedron* baseBoxPoly = aBox.GetPolyhedron();
120
121 if (query->queryMagnetic)
122 {
123 G4Colour boxZeroColourB = pointVisB->GetColor();
124 G4Colour boxZeroColourB2(boxZeroColourB.GetRed(), boxZeroColourB.GetGreen(), boxZeroColourB.GetBlue(), boxAlpha);
125 G4VisAttributes* boxZeroVisB = new G4VisAttributes(boxZeroColourB2);
126 G4String arrowPrefix = query->name + "_B_";
127 for (const auto& xyzBE: xyzBEs)
128 {
129 // x y z Bx By Bz Ex Ey Ez
130 G4ThreeVector B(xyzBE[3], xyzBE[4], xyzBE[5]);
131 G4ThreeVector unitB = B.unit();
132 G4double bMag = B.mag();
133 G4ThreeVector midPoint(xyzBE[0], xyzBE[1], xyzBE[2]);
134
135 // GetPolyhedron returns a pointer, but it's the same object always
136 // and there's no interface to allow rebuilding that polyhedron, so
137 // if we translate it we compound that transform. No method to just
138 // set the transform, so have to copy it.
139 G4Polyhedron boxPoly(*baseBoxPoly);
140 boxPoly.Transform(G4Transform3D(qRotation, midPoint));
141
142 if (bMag == 0)
143 {
144 sceneHandler.BeginPrimitives();
145 if (drawArrows)
146 {
147 G4Circle aPoint({xyzBE[0], xyzBE[1], xyzBE[2]});
148 aPoint.SetVisAttributes(pointVisB);
149 aPoint.SetDiameter(G4VMarker::SizeType::world, pointSize);
150 aPoint.SetFillStyle(G4VMarker::FillStyle::filled);
151 sceneHandler.AddPrimitive(aPoint);
152 }
153 if (drawBoxes)
154 {
155 boxPoly.SetVisAttributes(boxZeroVisB);
156 sceneHandler.AddPrimitive(boxPoly);
157 }
158 sceneHandler.EndPrimitives();
159 continue;
160 } // don't add 0 field arrows
161 G4ThreeVector startPoint = midPoint - 0.5*arrowLength*unitB;
162 G4ThreeVector endPoint = midPoint + 0.5*arrowLength*unitB;
163 G4double normalisedValue = bMag / maxFieldB;
164 if (!std::isfinite(normalisedValue))
165 {normalisedValue = 0.0;}
166 G4Colour arrowColour = bFieldColour.GetValue(normalisedValue);
167 if (drawBoxes)
168 {
169 G4Colour boxColour(arrowColour.GetRed(), arrowColour.GetGreen(), arrowColour.GetBlue(), boxAlpha);
170 sceneHandler.BeginPrimitives();
171 boxPoly.SetVisAttributes(G4VisAttributes(boxColour));
172 sceneHandler.AddPrimitive(boxPoly);
173 sceneHandler.EndPrimitives();
174 }
175 if (drawArrows)
176 {
177 std::stringstream ss;
178 ss << std::setw(2) << std::fixed << xyzBE[0] << "_" << xyzBE[1] << "_" << xyzBE[2];
179 std::string arrowName = ss.str();
180 G4ArrowModel FArrow(startPoint.x(), startPoint.y(), startPoint.z(),
181 endPoint.x(), endPoint.y(), endPoint.z(),
182 arrowWidth, arrowColour, arrowName);
183 FArrow.DescribeYourselfTo(sceneHandler);
184 }
185 }
186 delete boxZeroVisB;
187 }
188
189 if (query->queryElectric)
190 {
191 G4Colour boxZeroColourE = pointVisE->GetColor();
192 G4Colour boxZeroColourE2(boxZeroColourE.GetRed(),
193 boxZeroColourE.GetGreen(),
194 boxZeroColourE.GetBlue(),
195 boxAlpha);
196 G4VisAttributes* boxZeroVisE = new G4VisAttributes(boxZeroColourE2);
197 G4String arrowPrefix = query->name + "_E_";
198 for (const auto& xyzBE: xyzBEs)
199 {
200 // x y z Bx By Bz Ex Ey Ez
201 G4ThreeVector E(xyzBE[6], xyzBE[7], xyzBE[8]);
202 G4ThreeVector unitE = E.unit();
203 G4double eMag = E.mag();
204 G4ThreeVector midPoint(xyzBE[0], xyzBE[1], xyzBE[2]);
205
206 G4Polyhedron boxPoly(*baseBoxPoly);
207 boxPoly.Transform(G4Transform3D(qRotation, midPoint));
208 if (eMag == 0)
209 {
210 sceneHandler.BeginPrimitives();
211 if (drawArrows)
212 {
213 G4Circle aPoint({xyzBE[0], xyzBE[1], xyzBE[2]});
214 aPoint.SetVisAttributes(pointVisE);
215 aPoint.SetDiameter(G4VMarker::SizeType::world, pointSize);
216 aPoint.SetFillStyle(G4VMarker::FillStyle::filled);
217 sceneHandler.AddPrimitive(aPoint);
218 }
219 if (drawBoxes)
220 {
221 boxPoly.SetVisAttributes(boxZeroVisE);
222 sceneHandler.AddPrimitive(boxPoly);
223 }
224 sceneHandler.EndPrimitives();
225 continue;
226 } // don't add 0 field arrows
227 G4ThreeVector startPoint = midPoint - 0.5*arrowLength*unitE;
228 G4ThreeVector endPoint = midPoint + 0.5*arrowLength*unitE;
229 G4double normalisedValue = eMag / maxFieldE;
230 if (!std::isfinite(normalisedValue))
231 {normalisedValue = 0.0;}
232 G4Colour arrowColour = eFieldColour.GetValue(normalisedValue);
233 if (drawBoxes)
234 {
235 G4Colour boxColour(arrowColour.GetRed(), arrowColour.GetGreen(), arrowColour.GetBlue(), boxAlpha);
236 sceneHandler.BeginPrimitives();
237 boxPoly.SetVisAttributes(G4VisAttributes(boxColour));
238 sceneHandler.AddPrimitive(boxPoly);
239 sceneHandler.EndPrimitives();
240 }
241 if (drawArrows)
242 {
243 std::stringstream ss;
244 ss << std::setw(2) << std::fixed << xyzBE[0] << "_" << xyzBE[1] << "_" << xyzBE[2];
245 std::string arrowName = ss.str();
246 G4ArrowModel FArrow(startPoint.x(), startPoint.y(), startPoint.z(),
247 endPoint.x(), endPoint.y(), endPoint.z(),
248 arrowWidth, arrowColour, arrowName);
249 FArrow.DescribeYourselfTo(sceneHandler);
250 }
251 }
252 delete boxZeroVisE;
253 }
254
255 querier.CleanUp();
256 }
257}
258
259G4double BDSVisFieldModel::CalculateArrowLength(const BDSFieldQueryInfo* query) const
260{
261 G4double result = std::min({QIL(query->xInfo), QIL(query->yInfo), QIL(query->zInfo)});
262 return result;
263}
264
265G4double BDSVisFieldModel::QIL(const BDSFieldQueryInfo::QueryDimensionInfo& qi) const
266{
267 if (qi.n == 1)
268 {return std::numeric_limits<G4double>::max();}
269
270 G4double step = std::abs((qi.max - qi.min) / ((G4double)qi.n));
271 return 0.8 * step;
272}
273
274G4ThreeVector BDSVisFieldModel::BoxHalfSize(const BDSFieldQueryInfo* qi) const
275{
276 G4double minSize = std::numeric_limits<double>::max();
277 G4ThreeVector result(0,0,0);
278 std::array<G4bool,3> nonZero = {false, false, false};
279 std::array<const BDSFieldQueryInfo::QueryDimensionInfo,3> dims = {qi->xInfo, qi->yInfo, qi->zInfo};
280
281 G4int nDims = 0;
282 for (G4int i = 0; i < 3; i++)
283 {
284 nonZero[i] = dims[i].n > 1;
285 if (nonZero[i])
286 {
287 nDims += 1;
288 result[i] = dims[i].StepSize();
289 minSize = std::min(minSize, result[i]);
290 }
291 }
292
293 G4double thinSize;
294 if (nDims == 1)
295 {thinSize = 0.1 * minSize;} // 20% x 0.5
296 else // 2 or 3 dimensions - only poosibly set 3rd dimension as thin-ish
297 {thinSize = 1e-2 * minSize;}
298 for (G4int i = 0; i < 3; i++)
299 {result[i] = nonZero[i] ? result[i]*0.5 : thinSize;}
300
301 return result;
302}
BDSColourScaleMagma
Colour scale based on magma colour map.
Definition BDSColourScaleMagma.hh:49
BDSColourScaleViridis
Colour scale based on viridis colour map.
Definition BDSColourScaleViridis.hh:49
BDSFieldQueryForVis
Wrapper class of BDSFieldQuery that accumulates values for visualisation.
Definition BDSFieldQueryForVis.hh:37
BDSFieldQueryForVis::QueryField
virtual void QueryField(const BDSFieldQueryInfo *query)
Reserve the vector size we need first, then proceed as normal. Avoids big copies.
Definition BDSFieldQueryForVis.cc:48
BDSFieldQueryForVis::CleanUp
virtual void CleanUp()
Call the base class method and also clear the vector of values in this class.
Definition BDSFieldQueryForVis.cc:40
BDSFieldQueryInfo
Holder class for all information required for a field query.
Definition BDSFieldQueryInfo.hh:36
BDSVisFieldModel::CalculateArrowLength
G4double CalculateArrowLength(const BDSFieldQueryInfo *query) const
Return the minimum of 0.8 x the step length in x,y,z.
Definition BDSVisFieldModel.cc:259
BDSVisFieldModel::queries
const std::vector< BDSFieldQueryInfo * > queries
Cache of queries.
Definition BDSVisFieldModel.hh:63
BDSVisFieldModel::BoxHalfSize
G4ThreeVector BoxHalfSize(const BDSFieldQueryInfo *qi) const
Definition BDSVisFieldModel.cc:274
BDSVisFieldModel::QIL
G4double QIL(const BDSFieldQueryInfo::QueryDimensionInfo &qi) const
Query Info Length. For one dimension, return 0.8 x step size.
Definition BDSVisFieldModel.cc:265
BDSVisFieldModel::pointVisB
G4VisAttributes * pointVisB
The vis attributes for a zero-field marker.
Definition BDSVisFieldModel.hh:66
BDSVisFieldModel::pointVisE
G4VisAttributes * pointVisE
The vis attributes for a zero-field marker.
Definition BDSVisFieldModel.hh:67
BDSVisFieldModel::instanceCounter
static G4int instanceCounter
Definition BDSVisFieldModel.hh:48
Generated on Thu Nov 12 2026 21:37:00 for BDSIM by doxygen 1.9.8