Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 94206 2015-11-09 08:11: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 "G4Polyline.hh"
43 #include "G4VisAttributes.hh"
44 #include "G4SystemOfUnits.hh"
45 
46 #include <sstream>
47 #include <limits>
48 #include <vector>
49 
51 
53 (G4int nDataPointsPerMaxHalfScene,
54  Representation representation)
55 : fNDataPointsPerMaxHalfScene(nDataPointsPerMaxHalfScene)
56 , fRepresentation(representation)
57 {
58  fType = "G4MagneticFieldModel";
59  fGlobalTag = fType;
60  std::ostringstream oss;
61  oss << fNDataPointsPerMaxHalfScene;
62  fGlobalDescription = fType + ':' + oss.str();
63 }
64 
66 {
67 // G4cout << "G4MagneticFieldModel::DescribeYourselfTo" << G4endl;
68 
69  const G4VisExtent& extent = sceneHandler.GetExtent();
70  const G4double xMin = extent.GetXmin();
71  const G4double yMin = extent.GetYmin();
72  const G4double zMin = extent.GetZmin();
73  const G4double xMax = extent.GetXmax();
74  const G4double yMax = extent.GetYmax();
75  const G4double zMax = extent.GetZmax();
76  const G4double xHalfScene = 0.5 * (xMax - xMin);
77  const G4double yHalfScene = 0.5 * (yMax - yMin);
78  const G4double zHalfScene = 0.5 * (zMax - zMin);
79  const G4double xSceneCentre = 0.5 * (xMax + xMin);
80  const G4double ySceneCentre = 0.5 * (yMax + yMin);
81  const G4double zSceneCentre = 0.5 * (zMax + zMin);
82  const G4double maxHalfScene =
83  std::max(xHalfScene,std::max(yHalfScene,zHalfScene));
84  if (maxHalfScene <= 0.) {
85  G4cout
86  << "Extent non-positive."
87  << G4endl;
88  return;
89  }
90 
93  assert(tMgr);
95  assert(navigator);
96 
97  G4FieldManager* globalFieldMgr = tMgr->GetFieldManager();
98  const G4Field* globalField = 0;
99  const G4String intro = "G4MagneticFieldModel::DescribeYourselfTo: ";
100  if (globalFieldMgr) {
101  if (globalFieldMgr->DoesFieldExist()) {
102  globalField = globalFieldMgr->GetDetectorField();
103  if (!globalField) {
104  static G4bool warned = false;
105  if (!warned) {
106  G4cout << intro
107  << "Null global field pointer."
108  << G4endl;
109  warned = true;
110  }
111  }
112  } else {
113  static G4bool warned = false;
114  if (!warned) {
115  G4cout << intro
116  << "No global field exists."
117  << G4endl;
118  warned = true;
119  }
120  }
121  } else {
122  static G4bool warned = false;
123  if (!warned) {
124  G4cout << intro
125  << "No global field manager."
126  << G4endl;
127  warned = true;
128  }
129  }
130 
131  // Constants
132  const G4double interval = maxHalfScene / fNDataPointsPerMaxHalfScene;
133  const G4int nDataPointsPerXHalfScene = G4int(xHalfScene / interval);
134  const G4int nDataPointsPerYHalfScene = G4int(yHalfScene / interval);
135  const G4int nDataPointsPerZHalfScene = G4int(zHalfScene / interval);
136  const G4int nXSamples = 2 * nDataPointsPerXHalfScene + 1;
137  const G4int nYSamples = 2 * nDataPointsPerYHalfScene + 1;
138  const G4int nZSamples = 2 * nDataPointsPerZHalfScene + 1;
139  const G4int nSamples = nXSamples * nYSamples * nZSamples;
140  const G4int nSamples3 = nSamples * 3;
141  const G4double arrowLengthMax = 0.8 * interval;
142  const G4int nResults = 6; // 3 B-field + 3 E-field.
143 
144  // Working space for GetFieldValue.
145  double position_time[4] = {0,0,0,0};
146  double result[nResults];
147 
148  // Working vectors for field values, etc.
149  std::vector<G4double> BField(nSamples3); // Initialises to zero.
150  std::vector<G4double> BFieldMagnitude(nSamples); // Initialises to zero.
151  std::vector<G4double> xyz(nSamples3); // Initialises to zero.
152 
153  // Get field values and ascertain maximum field.
154  G4double BFieldMagnitudeMax = -std::numeric_limits<G4double>::max();
155  for (G4int i = 0; i < nXSamples; i++) {
156  G4double x = xSceneCentre + (i - nDataPointsPerXHalfScene) * interval;
157  position_time[0] = x;
158  for (G4int j = 0; j < nYSamples; j++) {
159  G4double y = ySceneCentre + (j - nDataPointsPerYHalfScene) * interval;
160  position_time[1] = y;
161  for (G4int k = 0; k < nZSamples; k++) {
162  G4double z = zSceneCentre + (k - nDataPointsPerZHalfScene) * interval;
163  position_time[2] = z;
164  // Calculate indices into working vectors
165  const G4int ijk = i * nYSamples * nZSamples + j * nZSamples + k;
166  const G4int ijk3 = ijk * 3;
167  // Find volume at this location.
168  G4ThreeVector pos(x,y,z);
169  const G4VPhysicalVolume* pPV =
170  navigator->LocateGlobalPointAndSetup(pos,0,false,true);
171  const G4Field* field = globalField;
172  if (pPV) {
173  // Get logical volume.
174  const G4LogicalVolume* pLV = pPV->GetLogicalVolume();
175  if (pLV) {
176  // Value for Region, if any, overrides
177  G4Region* pRegion = pLV->GetRegion();
178  if (pRegion) {
179  G4FieldManager* pRegionFieldMgr = pRegion->GetFieldManager();
180  if (pRegionFieldMgr) {
181  field = pRegionFieldMgr->GetDetectorField();
182  // G4cout << "Region with field" << G4endl;
183  }
184  }
185  // 'Local' value from logical volume, if any, overrides
186  G4FieldManager* pLVFieldMgr = pLV->GetFieldManager();
187  if (pLVFieldMgr) {
188  field = pLVFieldMgr->GetDetectorField();
189  // G4cout << "Logical volume with field" << G4endl;
190  }
191  }
192  }
193  // If field found, get values and store in working vectors.
194  if (field) {
195  // Get field values in result array.
196  field->GetFieldValue(position_time,result);
197  // G4cout
198  // << "BField/T:"
199  // << " " << result[0]/tesla
200  // << " " << result[1]/tesla
201  // << " " << result[2]/tesla
202  // << G4endl;
203  // Store B-field components.
204  for (G4int l = 0; l < 3; l++) {
205  BField[ijk3 + l] = result[l];
206  }
207  // Calculate magnitude and store.
208  G4double mag = sqrt
209  (result[0]*result[0]+result[1]*result[1]+result[2]*result[2]);
210  BFieldMagnitude[ijk] = mag;
211  // Store position.
212  xyz[ijk3] = x;
213  xyz[ijk3 + 1] = y;
214  xyz[ijk3 + 2] = z;
215  // Find maximum field magnitude.
216  if (mag > BFieldMagnitudeMax) {
217  BFieldMagnitudeMax = mag;
218  }
219  }
220  }
221  }
222  }
223 
224  if (BFieldMagnitudeMax <= 0) {
225  G4cout
226  << "No field in this scene."
227  << G4endl;
228  return;
229  }
230 
231  if (fRepresentation == Representation::lightArrow) sceneHandler.BeginPrimitives();
232  for (G4int i = 0; i < nSamples; i++) {
233  if (BFieldMagnitude[i] > 0) {
234  const G4int i3 = i * 3;
235  // Field (Bx,By,Bz) at (x,y,z).
236  const G4double Bx = BField[i3];
237  const G4double By = BField[i3 + 1];
238  const G4double Bz = BField[i3 + 2];
239  const G4double x = xyz[i3];
240  const G4double y = xyz[i3 + 1];
241  const G4double z = xyz[i3 + 2];
242  const G4double B = BFieldMagnitude[i];
243 // G4cout
244 // << "Position/mm, BField/T unpacked:"
245 // << ' ' << x/mm
246 // << ' ' << y/mm
247 // << ' ' << z/mm
248 // << " " << Bx/tesla
249 // << " " << By/tesla
250 // << " " << Bz/tesla
251 // << G4endl;
252  if (B > 0.) {
253  const G4double f = B / BFieldMagnitudeMax;
254  G4double red = 0., green = 0., blue = 0., alpha = 1.;
255  if (f < 0.5) { // Linear colour scale: 0->0.5->1 is blue->green->red.
256  green = 2. * f;
257  blue = 2. * (0.5 - f);
258  } else {
259  red = 2. * (f - 0.5);
260  green = 2. * (1.0 - f);
261  }
262  const G4Colour arrowColour(red,green,blue,alpha);
263  const G4double arrowLength = arrowLengthMax * f;
264  // Base of arrow is at (x,y,z).
265  const G4double& x1 = x;
266  const G4double& y1 = y;
267  const G4double& z1 = z;
268  // Head of arrow depends on field direction and strength.
269  const G4double x2 = x1 + arrowLength * Bx / B;
270  const G4double y2 = y1 + arrowLength * By / B;
271  const G4double z2 = z1 + arrowLength * Bz / B;
272  if (fRepresentation == Representation::fullArrow) {
273  G4ArrowModel BArrow(x1,y1,z1,x2,y2,z2,arrowLength/5,arrowColour);
274  BArrow.DescribeYourselfTo(sceneHandler);
275  } else if (fRepresentation == Representation::lightArrow) {
276  G4Polyline BArrowLite;
277  G4VisAttributes va(arrowColour);
278  BArrowLite.SetVisAttributes(va);
279  BArrowLite.push_back(G4Point3D(x1,y1,z1));
280  BArrowLite.push_back(G4Point3D(x2,y2,z2));
281  sceneHandler.AddPrimitive(BArrowLite);
282  }
283  }
284  }
285  }
286  if (fRepresentation == Representation::lightArrow) sceneHandler.EndPrimitives();
287 }
G4double G4ParticleHPJENDLHEData::G4double result
Definition: test07.cc:36
G4double GetXmin() const
Definition: G4VisExtent.hh:89
virtual void GetFieldValue(const double Point[4], double *fieldArr) const =0
G4double GetXmax() const
Definition: G4VisExtent.hh:90
G4Navigator * GetNavigatorForTracking() const
virtual void BeginPrimitives(const G4Transform3D &objectTransformation=G4Transform3D())=0
HepGeom::Point3D< G4double > G4Point3D
Definition: G4Point3D.hh:35
tuple x
Definition: test.py:50
G4Region * GetRegion() const
double B(double temperature)
Definition: test07.cc:36
int G4int
Definition: G4Types.hh:78
virtual void AddPrimitive(const G4Polyline &)=0
G4MagneticFieldModel(G4int nDataPointsPerHalfScene=10, Representation representation=Representation::fullArrow)
virtual void DescribeYourselfTo(G4VGraphicsScene &)
G4GLOB_DLL std::ostream G4cout
G4double GetYmax() const
Definition: G4VisExtent.hh:92
bool G4bool
Definition: G4Types.hh:79
G4bool DoesFieldExist() const
G4double GetZmax() const
Definition: G4VisExtent.hh:94
tuple navigator
Definition: write_gdml.py:35
void SetVisAttributes(const G4VisAttributes *)
Definition: G4Visible.cc:80
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:125
G4double GetZmin() const
Definition: G4VisExtent.hh:93
tuple z
Definition: test.py:28
G4double GetYmin() const
Definition: G4VisExtent.hh:91
#define G4endl
Definition: G4ios.hh:61
G4FieldManager * GetFieldManager() const
virtual void EndPrimitives()=0
virtual void DescribeYourselfTo(G4VGraphicsScene &)
double G4double
Definition: G4Types.hh:76
const G4Field * GetDetectorField() const
static const G4double alpha
static const G4double pos