ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VFieldModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VFieldModel.cc
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"
47 #include "G4VGraphicsScene.hh"
48 #include "G4VPhysicalVolume.hh"
49 #include "G4VisAttributes.hh"
50 
51 #include <sstream>
52 #include <limits>
53 #include <vector>
54 
55 
56 // Constructor and destructor
57 
59 
61 (const G4String& typeOfField, const G4String& symbol,
62  const G4VisExtent& extentForField,
63  const std::vector<G4PhysicalVolumesSearchScene::Findings>& pvFindings,
64  G4int nDataPointsPerMaxHalfScene,
65  Representation representation,
66  G4int arrow3DLineSegmentsPerCircle)
67 : fExtentForField(extentForField)
68 , fPVFindings(pvFindings)
69 , fNDataPointsPerMaxHalfScene(nDataPointsPerMaxHalfScene)
70 , fRepresentation(representation)
71 , fArrow3DLineSegmentsPerCircle(arrow3DLineSegmentsPerCircle)
72 , fTypeOfField(typeOfField)
73 , fArrowPrefix(symbol)
74 {
75  fType = "G4"+typeOfField+"FieldModel";
76  fGlobalTag = fType;
77 
78  std::ostringstream oss;
79  oss << ':' << fNDataPointsPerMaxHalfScene
80  << ':' << fArrow3DLineSegmentsPerCircle;
81  if (fExtentForField == G4VisExtent::GetNullExtent()) {
82  oss << " whole scene";
83  } else {
84  oss
85  << ':' << fExtentForField.GetXmin()
86  << ':' << fExtentForField.GetXmax()
87  << ':' << fExtentForField.GetYmin()
88  << ':' << fExtentForField.GetYmax()
89  << ':' << fExtentForField.GetZmin()
90  << ':' << fExtentForField.GetZmax();
91  }
92  for (const auto& findings: fPVFindings) {
93  oss
94  << ',' << findings.fpFoundPV->GetName()
95  << ':' << findings.fFoundPVCopyNo;
96  }
97  if (fRepresentation == Representation::fullArrow) {
98  oss << " full arrow";
99  } else if (fRepresentation == Representation::lightArrow) {
100  oss << " light arrow";
101  }
102 
103  fGlobalDescription = fType + oss.str();
104 }
105 
106 
107 // The main task of a model is to describe itself to the graphics scene.
108 
110 // G4cout << "G4VFieldModel::DescribeYourselfTo" << G4endl;
111 
114  assert(tMgr);
116  assert(navigator);
117 
118  G4FieldManager* globalFieldMgr = tMgr->GetFieldManager();
119  const G4Field* globalField = 0;
120  const G4String intro = "G4VFieldModel::DescribeYourselfTo: ";
121  if (globalFieldMgr) {
122  if (globalFieldMgr->DoesFieldExist()) {
123  globalField = globalFieldMgr->GetDetectorField();
124  if (!globalField) {
125  static G4bool warned = false;
126  if (!warned) {
127  G4cout << intro << "Null global field pointer." << G4endl;
128  warned = true;
129  }
130  }
131  }
132  } else {
133  static G4bool warned = false;
134  if (!warned) {
135  G4cout << intro << "No global field manager." << G4endl;
136  warned = true;
137  }
138  }
139 
140  G4VisExtent sceneExtent = sceneHandler.GetExtent();
141  const G4double& xMin = sceneExtent.GetXmin();
142  const G4double& yMin = sceneExtent.GetYmin();
143  const G4double& zMin = sceneExtent.GetZmin();
144  const G4double& xMax = sceneExtent.GetXmax();
145  const G4double& yMax = sceneExtent.GetYmax();
146  const G4double& zMax = sceneExtent.GetZmax();
147  const G4double xHalfScene = 0.5 * (xMax - xMin);
148  const G4double yHalfScene = 0.5 * (yMax - yMin);
149  const G4double zHalfScene = 0.5 * (zMax - zMin);
150  const G4double xSceneCentre = 0.5 * (xMax + xMin);
151  const G4double ySceneCentre = 0.5 * (yMax + yMin);
152  const G4double zSceneCentre = 0.5 * (zMax + zMin);
153  const G4double maxHalfScene =
154  std::max(xHalfScene,std::max(yHalfScene,zHalfScene));
155  if (maxHalfScene <= 0.) {
156  G4cout << "Scene extent non-positive." << G4endl;
157  return;
158  }
159 
160  // Constants
161  const G4double interval = maxHalfScene / fNDataPointsPerMaxHalfScene;
162  const G4int nDataPointsPerXHalfScene = G4int(xHalfScene / interval);
163  const G4int nDataPointsPerYHalfScene = G4int(yHalfScene / interval);
164  const G4int nDataPointsPerZHalfScene = G4int(zHalfScene / interval);
165  const G4int nXSamples = 2 * nDataPointsPerXHalfScene + 1;
166  const G4int nYSamples = 2 * nDataPointsPerYHalfScene + 1;
167  const G4int nZSamples = 2 * nDataPointsPerZHalfScene + 1;
168  const G4int nSamples = nXSamples * nYSamples * nZSamples;
169  const G4double arrowLengthMax = 0.8 * interval;
170 
171  // Working vectors for field values, etc.
172  std::vector<G4Point3D> Field(nSamples); // Initialises to (0,0,0)
173  std::vector<G4Point3D> xyz(nSamples); // Initialises to (0,0,0)
174  G4double FieldMagnitudeMax = -std::numeric_limits<G4double>::max();
175 
176  // Get field values and ascertain maximum field.
177  for (G4int i = 0; i < nXSamples; i++) {
178  G4double x = xSceneCentre + (i - nDataPointsPerXHalfScene) * interval;
179 
180  for (G4int j = 0; j < nYSamples; j++) {
181  G4double y = ySceneCentre + (j - nDataPointsPerYHalfScene) * interval;
182 
183  for (G4int k = 0; k < nZSamples; k++) {
184  G4double z = zSceneCentre + (k - nDataPointsPerZHalfScene) * interval;
185 
186  // Calculate indices into working vectors
187  const G4int ijk = i * nYSamples * nZSamples + j * nZSamples + k;
188  xyz[ijk].set(x,y,z);
189 
190  G4ThreeVector pos(x,y,z);
191 
192  // Check if point is in extent for field
194  const auto& ext = fExtentForField; // Alias
195  if (x < ext.GetXmin() || x > ext.GetXmax() ||
196  y < ext.GetYmin() || y > ext.GetYmax() ||
197  z < ext.GetZmin() || z > ext.GetZmax())
198  continue;
199  }
200 
201  // Check if point is in findings
202  if (!fPVFindings.empty()) {
203  G4bool isInPV = false;
204  for (const auto& findings: fPVFindings) {
205  G4VPhysicalVolume* pv = findings.fpFoundPV;
206  G4int copyNo = findings.fFoundPVCopyNo;
207  G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
208  G4PVParameterised* pvParam = dynamic_cast<G4PVParameterised*>(pv);
209  if (pvParam) {
210  auto* param = pvParam->GetParameterisation();
211  solid = param->ComputeSolid(copyNo,pvParam);
212  solid->ComputeDimensions(param,copyNo,pvParam);
213  }
214  // Transform point to local coordinate system
215  const auto& transform = findings.fFoundObjectTransformation;
216  auto rotation = transform.getRotation();
217  auto translation = transform.getTranslation();
218  G4ThreeVector lPos = pos; lPos -= translation; lPos.transform(rotation.invert());
219  if (solid->Inside(lPos)==kInside) {
220  isInPV = true;
221  break;
222  }
223  }
224  if (!isInPV) continue;
225  }
226  // Point is in findings - or there were no findings
227 
228  // Find volume and field at this location.
229  const G4VPhysicalVolume* pPV =
230  navigator->LocateGlobalPointAndSetup(pos,0,false,true);
231  const G4Field* field = globalField;
232  if (pPV) {
233  // Get logical volume.
234  const G4LogicalVolume* pLV = pPV->GetLogicalVolume();
235  if (pLV) {
236  // Value for Region, if any, overrides
237  G4Region* pRegion = pLV->GetRegion();
238  if (pRegion) {
239  G4FieldManager* pRegionFieldMgr = pRegion->GetFieldManager();
240  if (pRegionFieldMgr) {
241  field = pRegionFieldMgr->GetDetectorField();
242  // G4cout << "Region with field" << G4endl;
243  }
244  }
245  // 'Local' value from logical volume, if any, overrides
246  G4FieldManager* pLVFieldMgr = pLV->GetFieldManager();
247  if (pLVFieldMgr) {
248  field = pLVFieldMgr->GetDetectorField();
249  // G4cout << "Logical volume with field" << G4endl;
250  }
251  }
252  }
253 
254  G4double time = 0.; // FIXME: Can we get event time in some way?
255 
256  // Subclasses will have implemented this for their own field
257  GetFieldAtLocation(field, xyz[ijk], time, Field[ijk]);
258 
259  G4double mag = Field[ijk].mag();
260  if (mag > FieldMagnitudeMax) FieldMagnitudeMax = mag;
261  } // for (k, z
262  } // for (j, y
263  } // for (i, x
264 
265  if (FieldMagnitudeMax <= 0.) {
266  G4cout << "No " << fTypeOfField << " field in this extent." << G4endl;
267  return;
268  }
269 
270  for (G4int i = 0; i < nSamples; i++) {
271  const G4double Fmag = Field[i].mag();
272  const G4double f = Fmag / FieldMagnitudeMax;
273  if (f <= 0.) continue; // Skip zero field locations
274 
275  G4double red = 0., green = 0., blue = 0., alpha = 1.;
276  if (f < 0.5) { // Linear colour scale: 0->0.5->1 is red->green->blue.
277  green = 2. * f;
278  red = 2. * (0.5 - f);
279  } else {
280  blue = 2. * (f - 0.5);
281  green = 2. * (1.0 - f);
282  }
283  const G4Colour arrowColour(red,green,blue,alpha);
284 
285  // Very small arrows are difficult to see. Better to draw a line.
286  G4bool drawAsLine = false;
287  switch (fRepresentation) {
288  case Representation::fullArrow:
289  if (f < 0.1) {
290  drawAsLine = true;
291  }
292  break;
293  case Representation::lightArrow:
294  drawAsLine = true;
295  break;
296  default:
297  break;
298  }
299 
300  // Head of arrow depends on field direction and strength...
301  G4double arrowLength = arrowLengthMax * f;
302  // ...but limit the length so it's visible.
303  if (f < 0.01) arrowLength = arrowLengthMax * 0.01;
304  const G4Point3D head = xyz[i] + arrowLength*Field[i]/Fmag;
305 
306  if (drawAsLine) {
307  G4Polyline FArrowLite;
308  G4VisAttributes va(arrowColour);
309  va.SetLineWidth(2.);
310  FArrowLite.SetVisAttributes(va);
311  FArrowLite.push_back(xyz[i]);
312  FArrowLite.push_back(head);
313  sceneHandler.BeginPrimitives();
314  sceneHandler.AddPrimitive(FArrowLite);
315  sceneHandler.EndPrimitives();
316  } else {
317  G4ArrowModel FArrow(xyz[i].x(), xyz[i].y(), xyz[i].z(),
318  head.x(), head.y(), head.z(),
319  arrowLength/5, arrowColour,
320  fArrowPrefix+"Field",
322  FArrow.DescribeYourselfTo(sceneHandler);
323  }
324  } // for (i, nSamples
325 }