| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Michael Kelsey 31 January 2019
27 //
28 // Class Description:
29 //
30 // Abstract base class to implement drawing vector field geometries
31 // (e.g., electric, magnetic or gravity). Implementation extracted
32 // from G4MagneticFieldModel, with field-value access left pure
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& symbol,
63 const G4VisExtent& extentForField,
64 const std::vector<G4PhysicalVolumesSearchScene::Findings>& pvFindings,
65 G4int nDataPointsPerMaxHalfExtent,
66 Representation representation,
67 G4int arrow3DLineSegmentsPerCircle)
68 : fExtentForField(extentForField)
69 , fPVFindings(pvFindings)
70 , fNDataPointsPerMaxHalfExtent(nDataPointsPerMaxHalfExtent)
71 , fRepresentation(representation)
72 , fArrow3DLineSegmentsPerCircle(arrow3DLineSegmentsPerCircle)
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::GetNullExtent()) {
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::fullArrow) {
99 oss << " full arrow";
100 } else if (fRepresentation == Representation::lightArrow) {
101 oss << " light arrow";
102 }
103
104 fGlobalDescription = fType + oss.str();
105 }
106
107
108 // The main task of a model is to describe itself to the graphics scene.
109
110 void G4VFieldModel::DescribeYourselfTo(G4VGraphicsScene& sceneHandler) {
111 // G4cout << "G4VFieldModel::DescribeYourselfTo" << G4endl;
112
113 G4TransportationManager* tMgr =
114 G4TransportationManager::GetTransportationManager();
115 assert(tMgr);
116 G4Navigator* navigator = tMgr->GetNavigatorForTracking();
117 assert(navigator);
118
119 G4FieldManager* globalFieldMgr = tMgr->GetFieldManager();
120 const G4Field* globalField = 0;
121 const G4String intro = "G4VFieldModel::DescribeYourselfTo: ";
122 if (globalFieldMgr) {
123 if (globalFieldMgr->DoesFieldExist()) {
124 globalField = globalFieldMgr->GetDetectorField();
125 if (!globalField) {
126 static G4bool warned = false;
127 if (!warned) {
128 G4warn << intro << "Null global field pointer." << G4endl;
129 warned = true;
130 }
131 }
132 }
133 } else {
134 static G4bool warned = false;
135 if (!warned) {
136 G4warn << intro << "No global field manager." << G4endl;
137 warned = true;
138 }
139 }
140
141 G4VisExtent extent = sceneHandler.GetExtent();
142 if (fExtentForField == G4VisExtent::GetNullExtent()) {
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 - xMin);
154 const G4double yHalfScene = 0.5 * (yMax - yMin);
155 const G4double zHalfScene = 0.5 * (zMax - zMin);
156 const G4double xSceneCentre = 0.5 * (xMax + xMin);
157 const G4double ySceneCentre = 0.5 * (yMax + yMin);
158 const G4double zSceneCentre = 0.5 * (zMax + zMin);
159 const G4double maxHalfScene =
160 std::max(xHalfScene,std::max(yHalfScene,zHalfScene));
161 if (maxHalfScene <= 0.) {
162 G4warn << "Scene extent non-positive." << G4endl;
163 return;
164 }
165
166 // Constants
167 const G4double interval = maxHalfScene / fNDataPointsPerMaxHalfExtent;
168 const G4int nDataPointsPerXHalfScene = G4int(xHalfScene / interval);
169 const G4int nDataPointsPerYHalfScene = G4int(yHalfScene / interval);
170 const G4int nDataPointsPerZHalfScene = G4int(zHalfScene / interval);
171 const G4int nXSamples = 2 * nDataPointsPerXHalfScene + 1;
172 const G4int nYSamples = 2 * nDataPointsPerYHalfScene + 1;
173 const G4int nZSamples = 2 * nDataPointsPerZHalfScene + 1;
174 const G4int nSamples = nXSamples * nYSamples * nZSamples;
175 const G4double arrowLengthMax = 0.8 * interval;
176
177 // Working vectors for field values, etc.
178 std::vector<G4Point3D> Field(nSamples); // Initialises to (0,0,0)
179 std::vector<G4Point3D> xyz(nSamples); // Initialises to (0,0,0)
180 G4double FieldMagnitudeMax = -std::numeric_limits<G4double>::max();
181
182 // Get field values and ascertain maximum field.
183 for (G4int i = 0; i < nXSamples; i++) {
184 G4double x = xSceneCentre + (i - nDataPointsPerXHalfScene) * interval;
185
186 for (G4int j = 0; j < nYSamples; j++) {
187 G4double y = ySceneCentre + (j - nDataPointsPerYHalfScene) * interval;
188
189 for (G4int k = 0; k < nZSamples; k++) {
190 G4double z = zSceneCentre + (k - nDataPointsPerZHalfScene) * interval;
191
192 // Calculate indices into working vectors
193 const G4int ijk = i * nYSamples * nZSamples + j * nZSamples + k;
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: fPVFindings) {
202 G4VPhysicalVolume* pv = findings.fpFoundPV;
203 G4int copyNo = findings.fFoundPVCopyNo;
204 G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
205 G4PVParameterised* pvParam = dynamic_cast<G4PVParameterised*>(pv);
206 if (pvParam) {
207 auto* param = pvParam->GetParameterisation();
208 solid = param->ComputeSolid(copyNo,pvParam);
209 solid->ComputeDimensions(param,copyNo,pvParam);
210 }
211 // Transform point to local coordinate system
212 const auto& transform = findings.fFoundObjectTransformation;
213 auto rotation = transform.getRotation();
214 auto translation = transform.getTranslation();
215 G4ThreeVector lPos = pos; lPos -= translation; lPos.transform(rotation.invert());
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 were no findings
224
225 // Find volume and field at this location.
226 const G4VPhysicalVolume* pPV =
227 navigator->LocateGlobalPointAndSetup(pos,0,false,true);
228 const G4Field* field = globalField;
229 if (pPV) {
230 // Get logical volume.
231 const G4LogicalVolume* pLV = pPV->GetLogicalVolume();
232 if (pLV) {
233 // Value for Region, if any, overrides
234 G4Region* pRegion = pLV->GetRegion();
235 if (pRegion) {
236 G4FieldManager* pRegionFieldMgr = pRegion->GetFieldManager();
237 if (pRegionFieldMgr) {
238 field = pRegionFieldMgr->GetDetectorField();
239 // G4cout << "Region with field" << G4endl;
240 }
241 }
242 // 'Local' value from logical volume, if any, overrides
243 G4FieldManager* pLVFieldMgr = pLV->GetFieldManager();
244 if (pLVFieldMgr) {
245 field = pLVFieldMgr->GetDetectorField();
246 // G4cout << "Logical volume with field" << G4endl;
247 }
248 }
249 }
250
251 G4double time = 0.; // FIXME: Can we get event time in some way?
252
253 // Subclasses will have implemented this for their own field
254 GetFieldAtLocation(field, xyz[ijk], time, Field[ijk]);
255
256 G4double mag = Field[ijk].mag();
257 if (mag > FieldMagnitudeMax) FieldMagnitudeMax = mag;
258 } // for (k, z
259 } // for (j, y
260 } // for (i, x
261
262 if (FieldMagnitudeMax <= 0.) {
263 G4warn << "No " << fTypeOfField << " field in this extent." << G4endl;
264 return;
265 }
266
267 for (G4int i = 0; i < nSamples; i++) {
268 const G4double Fmag = Field[i].mag();
269 const G4double f = Fmag / FieldMagnitudeMax;
270 if (f <= 0.) continue; // Skip zero field locations
271
272 G4double red = 0., green = 0., blue = 0., alpha = 1.;
273 if (f < 0.5) { // Linear colour scale: 0->0.5->1 is red->green->blue.
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,alpha);
281
282 // Very small arrows are difficult to see. Better to draw a line.
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 direction and strength...
298 G4double arrowLength = arrowLengthMax * f;
299 // ...but limit the length so it's visible.
300 if (f < 0.01) arrowLength = arrowLengthMax * 0.01;
301 const G4Point3D head = xyz[i] + arrowLength*Field[i]/Fmag;
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(), xyz[i].z(),
315 head.x(), head.y(), head.z(),
316 arrowLength/5, arrowColour,
317 fArrowPrefix+"Field",
318 fArrow3DLineSegmentsPerCircle);
319 FArrow.DescribeYourselfTo(sceneHandler);
320 }
321 } // for (i, nSamples
322 }
323