Geant4  10.00.p03
G4MagneticFieldModel.cc
Go to the documentation of this file.
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 //
27 // $Id: G4MagneticFieldModel.cc 74097 2013-09-22 16:03:59Z gcosmo $
28 //
29 //
30 // John Allison 17th August 2013
31 // Model that knows how to draw the magnetic field.
32 
33 #include "G4MagneticFieldModel.hh"
34 
35 #include "G4VGraphicsScene.hh"
37 #include "G4FieldManager.hh"
38 #include "G4Field.hh"
39 #include "G4Colour.hh"
40 #include "G4VPhysicalVolume.hh"
41 #include "G4ArrowModel.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include <sstream>
45 #include <limits>
46 #include <vector>
47 
49 {
50 }
51 
53 : fNDataPointsPerMaxHalfScene(nDataPointsPerMaxHalfScene)
54 {
55  fType = "G4MagneticFieldModel";
56  fGlobalTag = fType;
57  std::ostringstream oss;
59  fGlobalDescription = fType + ':' + oss.str();
60 }
61 
63 {
64 // G4cout << "G4MagneticFieldModel::DescribeYourselfTo" << G4endl;
65 
66  const G4VisExtent& extent = sceneHandler.GetExtent();
67  const G4double xMin = extent.GetXmin();
68  const G4double yMin = extent.GetYmin();
69  const G4double zMin = extent.GetZmin();
70  const G4double xMax = extent.GetXmax();
71  const G4double yMax = extent.GetYmax();
72  const G4double zMax = extent.GetZmax();
73  const G4double xHalfScene = 0.5 * (xMax - xMin);
74  const G4double yHalfScene = 0.5 * (yMax - yMin);
75  const G4double zHalfScene = 0.5 * (zMax - zMin);
76 // const G4double xSceneCentre = 0.5 * (xMax + xMin);
77 // const G4double ySceneCentre = 0.5 * (yMax + yMin);
78 // const G4double zSceneCentre = 0.5 * (zMax + zMin);
79  const G4double maxHalfScene =
80  std::max(xHalfScene,std::max(yHalfScene,zHalfScene));
81  if (maxHalfScene <= 0.) {
82  G4cout
83  << "Extent non-positive."
84  << G4endl;
85  return;
86  }
87 
90  assert(tMgr);
91  G4Navigator* navigator = tMgr->GetNavigatorForTracking();
92  assert(navigator);
93 
94  G4FieldManager* globalFieldMgr = tMgr->GetFieldManager();
95  const G4Field* globalField = 0;
96  if (globalFieldMgr) {
97  if (globalFieldMgr->DoesFieldExist()) {
98  globalField = globalFieldMgr->GetDetectorField();
99  if (!globalField) {
100  G4cout
101  << "Null global field pointer."
102  << G4endl;
103  }
104  } else {
105  G4cout
106  << "No global field exists."
107  << G4endl;
108  }
109  } else {
110  G4cout
111  << "No global field manager."
112  << G4endl;
113  }
114 
115  // Constants
116  const G4double interval = maxHalfScene / fNDataPointsPerMaxHalfScene;
117  const G4int nDataPointsPerXHalfScene = G4int(xHalfScene / interval);
118  const G4int nDataPointsPerYHalfScene = G4int(yHalfScene / interval);
119  const G4int nDataPointsPerZHalfScene = G4int(zHalfScene / interval);
120  const G4int nXSamples = 2 * nDataPointsPerXHalfScene + 1;
121  const G4int nYSamples = 2 * nDataPointsPerYHalfScene + 1;
122  const G4int nZSamples = 2 * nDataPointsPerZHalfScene + 1;
123  const G4int nSamples = nXSamples * nYSamples * nZSamples;
124  const G4int nSamples3 = nSamples * 3;
125  const G4double arrowLengthMax = 0.8 * interval;
126  const G4int nResults = 6; // 3 B-field + 3 E-field.
127 
128  // Working space for GetFieldValue.
129  double position_time[4] = {0,0,0,0};
130  double result[nResults];
131 
132  // Working vectors for field values, etc.
133  std::vector<G4double> BField(nSamples3); // Initialises to zero.
134  std::vector<G4double> BFieldMagnitude(nSamples); // Initialises to zero.
135  std::vector<G4double> xyz(nSamples3); // Initialises to zero.
136 
137  // Get field values and ascertain maximum field.
138  G4double BFieldMagnitudeMax = -std::numeric_limits<G4double>::max();
139  for (G4int i = 0; i < nXSamples; i++) {
140  G4double x = (i - nDataPointsPerXHalfScene) * interval;
141  position_time[0] = x;
142  for (G4int j = 0; j < nYSamples; j++) {
143  G4double y = (j - nDataPointsPerYHalfScene) * interval;
144  position_time[1] = y;
145  for (G4int k = 0; k < nZSamples; k++) {
146  G4double z = (k - nDataPointsPerZHalfScene) * interval;
147  position_time[2] = z;
148  // Calculate indices into working vectors
149  const G4int ijk = i * nYSamples * nZSamples + j * nZSamples + k;
150  const G4int ijk3 = ijk * 3;
151  // Find volume at this location.
152  G4ThreeVector pos(x,y,z);
153  const G4VPhysicalVolume* pPV =
154  navigator->LocateGlobalPointAndSetup(pos,0,false,true);
155  const G4Field* field = globalField;
156  if (pPV) {
157  // Get logical volume.
158  const G4LogicalVolume* pLV = pPV->GetLogicalVolume();
159  if (pLV) {
160  // Value for Region, if any, overrides
161  G4Region* pRegion = pLV->GetRegion();
162  if (pRegion) {
163  G4FieldManager* pRegionFieldMgr = pRegion->GetFieldManager();
164  if (pRegionFieldMgr) {
165  field = pRegionFieldMgr->GetDetectorField();
166  // G4cout << "Region with field" << G4endl;
167  }
168  }
169  // 'Local' value from logical volume, if any, overrides
170  G4FieldManager* pLVFieldMgr = pLV->GetFieldManager();
171  if (pLVFieldMgr) {
172  field = pLVFieldMgr->GetDetectorField();
173  // G4cout << "Logical volume with field" << G4endl;
174  }
175  }
176  }
177  // If field found, get values and store in working vectors.
178  if (field) {
179  // Get field values in result array.
180  field->GetFieldValue(position_time,result);
181  // G4cout
182  // << "BField/T:"
183  // << " " << result[0]/tesla
184  // << " " << result[1]/tesla
185  // << " " << result[2]/tesla
186  // << G4endl;
187  // Store B-field components.
188  for (G4int l = 0; l < 3; l++) {
189  BField[ijk3 + l] = result[l];
190  }
191  // Calculate magnitude and store.
192  G4double mag = sqrt
193  (result[0]*result[0]+result[1]*result[1]+result[2]*result[2]);
194  BFieldMagnitude[ijk] = mag;
195  // Store position.
196  xyz[ijk3] = x;
197  xyz[ijk3 + 1] = y;
198  xyz[ijk3 + 2] = z;
199  // Find maximum field magnitude.
200  if (mag > BFieldMagnitudeMax) {
201  BFieldMagnitudeMax = mag;
202  }
203  }
204  }
205  }
206  }
207 
208  if (BFieldMagnitudeMax <= 0) {
209  G4cout
210  << "No field in this scene."
211  << G4endl;
212  return;
213  }
214 
215  for (G4int i = 0; i < nSamples; i++) {
216  if (BFieldMagnitude[i] > 0) {
217  const G4int i3 = i * 3;
218  const G4double x = xyz[i3];
219  const G4double y = xyz[i3 + 1];
220  const G4double z = xyz[i3 + 2];
221  const G4double B = BFieldMagnitude[i];
222 // G4cout
223 // << "Position/mm, BField/T unpacked:"
224 // << ' ' << x/mm
225 // << ' ' << y/mm
226 // << ' ' << z/mm
227 // << " " << BField[i3]/tesla
228 // << " " << BField[i3 + 1]/tesla
229 // << " " << BField[i3 + 2]/tesla
230 // << G4endl;
231  const G4double f = B / BFieldMagnitudeMax;
232  const G4double arrowLength = arrowLengthMax * f;
233  G4double red = 0., green = 0., blue = 0., alpha = 1.;
234  if (f < 0.5) { // Linear colour scale: 0->0.5->1 is blue->green->red.
235  green = 2. * f;
236  blue = 2. * (0.5 - f);
237  } else {
238  red = 2. * (f - 0.5);
239  green = 2. * (1.0 - f);
240  }
241  const G4Colour arrowColour(red,green,blue,alpha);
242  G4ArrowModel BArrow
243  (x,y,z,
244  x + arrowLength * BField[i3] / B,
245  y + arrowLength * BField[i3 + 1] / B,
246  z + arrowLength * BField[i3 + 2] / B,
247  arrowLength/5,
248  arrowColour);
249  BArrow.DescribeYourselfTo(sceneHandler);
250  }
251  }
252 }
Definition: test07.cc:36
G4double GetXmin() const
Definition: G4VisExtent.hh:89
CLHEP::Hep3Vector G4ThreeVector
G4String fType
Definition: G4VModel.hh:108
virtual void GetFieldValue(const double Point[4], double *fieldArr) const =0
G4double z
Definition: TRTMaterials.hh:39
#define assert(x)
Definition: mymalloc.cc:1309
G4double GetXmax() const
Definition: G4VisExtent.hh:90
G4Navigator * GetNavigatorForTracking() const
G4Region * GetRegion() const
Definition: test07.cc:36
int G4int
Definition: G4Types.hh:78
G4String fGlobalTag
Definition: G4VModel.hh:109
G4FieldManager * GetFieldManager() const
virtual void DescribeYourselfTo(G4VGraphicsScene &)
G4GLOB_DLL std::ostream G4cout
G4double GetYmax() const
Definition: G4VisExtent.hh:92
G4bool DoesFieldExist() const
G4String fGlobalDescription
Definition: G4VModel.hh:110
G4double GetZmax() const
Definition: G4VisExtent.hh:94
static G4TransportationManager * GetTransportationManager()
G4FieldManager * GetFieldManager() const
G4FieldManager * GetFieldManager() const
virtual const G4VisExtent & GetExtent() const
G4LogicalVolume * GetLogicalVolume() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:118
G4double GetZmin() const
Definition: G4VisExtent.hh:93
G4double GetYmin() const
Definition: G4VisExtent.hh:91
#define G4endl
Definition: G4ios.hh:61
G4MagneticFieldModel(G4int nDataPointsPerHalfScene=10)
virtual void DescribeYourselfTo(G4VGraphicsScene &)
double G4double
Definition: G4Types.hh:76
const G4Field * GetDetectorField() const
static const G4double alpha
static const G4double pos