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