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