| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Michael Kelsey 31 January 2019
27 //
28 // Class Description:
29 //
30 // Abstract base class to implement drawing ve
31 // (e.g., electric, magnetic or gravity). Imp
32 // from G4MagneticFieldModel, with field-value
33 // virtual for implementation by base classes.
34
35 #include "G4VFieldModel.hh"
36
37 #include "G4ArrowModel.hh"
38 #include "G4Colour.hh"
39 #include "G4Field.hh"
40 #include "G4FieldManager.hh"
41 #include "G4PVPlacement.hh"
42 #include "G4PVParameterised.hh"
43 #include "G4Point3D.hh"
44 #include "G4Polyline.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4TransportationManager.hh"
47 #include "G4VGraphicsScene.hh"
48 #include "G4VPhysicalVolume.hh"
49 #include "G4VisAttributes.hh"
50
51 #include <sstream>
52 #include <limits>
53 #include <vector>
54
55 #define G4warn G4cout
56
57 // Constructor and destructor
58
59 G4VFieldModel::~G4VFieldModel() {;}
60
61 G4VFieldModel::G4VFieldModel
62 (const G4String& typeOfField, const G4String&
63 const G4VisExtent& extentForField,
64 const std::vector<G4PhysicalVolumesSearchScen
65 G4int nDataPointsPerMaxHalfExtent,
66 Representation representation,
67 G4int arrow3DLineSegmentsPerCircle)
68 : fExtentForField(extentForField)
69 , fPVFindings(pvFindings)
70 , fNDataPointsPerMaxHalfExtent(nDataPointsPerM
71 , fRepresentation(representation)
72 , fArrow3DLineSegmentsPerCircle(arrow3DLineSeg
73 , fTypeOfField(typeOfField)
74 , fArrowPrefix(symbol)
75 {
76 fType = "G4"+typeOfField+"FieldModel";
77 fGlobalTag = fType;
78
79 std::ostringstream oss;
80 oss << ':' << fNDataPointsPerMaxHalfExtent
81 << ':' << fArrow3DLineSegmentsPerCircle;
82 if (fExtentForField == G4VisExtent::GetNullE
83 oss << " whole scene";
84 } else {
85 oss
86 << ':' << fExtentForField.GetXmin()
87 << ':' << fExtentForField.GetXmax()
88 << ':' << fExtentForField.GetYmin()
89 << ':' << fExtentForField.GetYmax()
90 << ':' << fExtentForField.GetZmin()
91 << ':' << fExtentForField.GetZmax();
92 }
93 for (const auto& findings: fPVFindings) {
94 oss
95 << ',' << findings.fpFoundPV->GetName()
96 << ':' << findings.fFoundPVCopyNo;
97 }
98 if (fRepresentation == Representation::fullA
99 oss << " full arrow";
100 } else if (fRepresentation == Representation
101 oss << " light arrow";
102 }
103
104 fGlobalDescription = fType + oss.str();
105 }
106
107
108 // The main task of a model is to describe its
109
110 void G4VFieldModel::DescribeYourselfTo(G4VGrap
111 // G4cout << "G4VFieldModel::DescribeYourself
112
113 G4TransportationManager* tMgr =
114 G4TransportationManager::GetTransportationMa
115 assert(tMgr);
116 G4Navigator* navigator = tMgr->GetNavigatorF
117 assert(navigator);
118
119 G4FieldManager* globalFieldMgr = tMgr->GetFi
120 const G4Field* globalField = 0;
121 const G4String intro = "G4VFieldModel::Descr
122 if (globalFieldMgr) {
123 if (globalFieldMgr->DoesFieldExist()) {
124 globalField = globalFieldMgr->GetDetecto
125 if (!globalField) {
126 static G4bool warned = false;
127 if (!warned) {
128 G4warn << intro << "Null global fiel
129 warned = true;
130 }
131 }
132 }
133 } else {
134 static G4bool warned = false;
135 if (!warned) {
136 G4warn << intro << "No global field mana
137 warned = true;
138 }
139 }
140
141 G4VisExtent extent = sceneHandler.GetExtent(
142 if (fExtentForField == G4VisExtent::GetNullE
143 extent = sceneHandler.GetExtent();
144 } else {
145 extent = fExtentForField;
146 }
147 const G4double& xMin = extent.GetXmin();
148 const G4double& yMin = extent.GetYmin();
149 const G4double& zMin = extent.GetZmin();
150 const G4double& xMax = extent.GetXmax();
151 const G4double& yMax = extent.GetYmax();
152 const G4double& zMax = extent.GetZmax();
153 const G4double xHalfScene = 0.5 * (xMax - xM
154 const G4double yHalfScene = 0.5 * (yMax - yM
155 const G4double zHalfScene = 0.5 * (zMax - zM
156 const G4double xSceneCentre = 0.5 * (xMax +
157 const G4double ySceneCentre = 0.5 * (yMax +
158 const G4double zSceneCentre = 0.5 * (zMax +
159 const G4double maxHalfScene =
160 std::max(xHalfScene,std::max(yHalfScene,zHal
161 if (maxHalfScene <= 0.) {
162 G4warn << "Scene extent non-positive." <<
163 return;
164 }
165
166 // Constants
167 const G4double interval = maxHalfScene / fND
168 const G4int nDataPointsPerXHalfScene = G4int
169 const G4int nDataPointsPerYHalfScene = G4int
170 const G4int nDataPointsPerZHalfScene = G4int
171 const G4int nXSamples = 2 * nDataPointsPerXH
172 const G4int nYSamples = 2 * nDataPointsPerYH
173 const G4int nZSamples = 2 * nDataPointsPerZH
174 const G4int nSamples = nXSamples * nYSamples
175 const G4double arrowLengthMax = 0.8 * interv
176
177 // Working vectors for field values, etc.
178 std::vector<G4Point3D> Field(nSamples);
179 std::vector<G4Point3D> xyz(nSamples);
180 G4double FieldMagnitudeMax = -std::numeric_l
181
182 // Get field values and ascertain maximum fi
183 for (G4int i = 0; i < nXSamples; i++) {
184 G4double x = xSceneCentre + (i - nDataPoin
185
186 for (G4int j = 0; j < nYSamples; j++) {
187 G4double y = ySceneCentre + (j - nDataPo
188
189 for (G4int k = 0; k < nZSamples; k++) {
190 G4double z = zSceneCentre + (k - nData
191
192 // Calculate indices into working vect
193 const G4int ijk = i * nYSamples * nZSa
194 xyz[ijk].set(x,y,z);
195
196 G4ThreeVector pos(x,y,z);
197
198 // Check if point is in findings
199 if (!fPVFindings.empty()) {
200 G4bool isInPV = false;
201 for (const auto& findings: fPVFindin
202 G4VPhysicalVolume* pv = findings.f
203 G4int copyNo = findings.fFoundPVCo
204 G4VSolid* solid = pv->GetLogicalVo
205 G4PVParameterised* pvParam = dynam
206 if (pvParam) {
207 auto* param = pvParam->GetParame
208 solid = param->ComputeSolid(copy
209 solid->ComputeDimensions(param,c
210 }
211 // Transform point to local coordi
212 const auto& transform = findings.f
213 auto rotation = transform.getRotat
214 auto translation = transform.getTr
215 G4ThreeVector lPos = pos; lPos -=
216 if (solid->Inside(lPos)==kInside)
217 isInPV = true;
218 break;
219 }
220 }
221 if (!isInPV) continue;
222 }
223 // Point is in findings - or there wer
224
225 // Find volume and field at this locat
226 const G4VPhysicalVolume* pPV =
227 navigator->LocateGlobalPointAndSetup(p
228 const G4Field* field = globalField;
229 if (pPV) {
230 // Get logical volume.
231 const G4LogicalVolume* pLV = pPV->Ge
232 if (pLV) {
233 // Value for Region, if any, overr
234 G4Region* pRegion = pLV->GetRegion
235 if (pRegion) {
236 G4FieldManager* pRegionFieldMgr
237 if (pRegionFieldMgr) {
238 field = pRegionFieldMgr->GetDe
239 // G4cout << "Region with fiel
240 }
241 }
242 // 'Local' value from logical volu
243 G4FieldManager* pLVFieldMgr = pLV-
244 if (pLVFieldMgr) {
245 field = pLVFieldMgr->GetDetector
246 // G4cout << "Logical volume wit
247 }
248 }
249 }
250
251 G4double time = 0.; // FIXME: Can we get ev
252
253 // Subclasses will have implemented this for
254 GetFieldAtLocation(field, xyz[ijk], time, Fi
255
256 G4double mag = Field[ijk].mag();
257 if (mag > FieldMagnitudeMax) FieldMagnitudeM
258 } // for (k, z
259 } // for (j, y
260 } // for (i, x
261
262 if (FieldMagnitudeMax <= 0.) {
263 G4warn << "No " << fTypeOfField << " field
264 return;
265 }
266
267 for (G4int i = 0; i < nSamples; i++) {
268 const G4double Fmag = Field[i].mag();
269 const G4double f = Fmag / FieldMagnitudeMa
270 if (f <= 0.) continue; // Skip zero field
271
272 G4double red = 0., green = 0., blue = 0.,
273 if (f < 0.5) { // Linear colour scale: 0-
274 green = 2. * f;
275 red = 2. * (0.5 - f);
276 } else {
277 blue = 2. * (f - 0.5);
278 green = 2. * (1.0 - f);
279 }
280 const G4Colour arrowColour(red,green,blue,
281
282 // Very small arrows are difficult to see.
283 G4bool drawAsLine = false;
284 switch (fRepresentation) {
285 case Representation::fullArrow:
286 if (f < 0.1) {
287 drawAsLine = true;
288 }
289 break;
290 case Representation::lightArrow:
291 drawAsLine = true;
292 break;
293 default:
294 break;
295 }
296
297 // Head of arrow depends on field directio
298 G4double arrowLength = arrowLengthMax * f;
299 // ...but limit the length so it's visible
300 if (f < 0.01) arrowLength = arrowLengthMax
301 const G4Point3D head = xyz[i] + arrowLengt
302
303 if (drawAsLine) {
304 G4Polyline FArrowLite;
305 G4VisAttributes va(arrowColour);
306 va.SetLineWidth(2.);
307 FArrowLite.SetVisAttributes(va);
308 FArrowLite.push_back(xyz[i]);
309 FArrowLite.push_back(head);
310 sceneHandler.BeginPrimitives();
311 sceneHandler.AddPrimitive(FArrowLite);
312 sceneHandler.EndPrimitives();
313 } else {
314 G4ArrowModel FArrow(xyz[i].x(), xyz[i].y
315 head.x(), head.y(),
316 arrowLength/5, arrow
317 fArrowPrefix+"Field"
318 fArrow3DLineSegments
319 FArrow.DescribeYourselfTo(sceneHandler);
320 }
321 } // for (i, nSamples
322 }
323