| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4MagneticFieldModel.cc 94206 2015-11-09 08:11:59Z gcosmo $
27 // 28 //
28 // 29 //
29 // John Allison 17th August 2013 30 // John Allison 17th August 2013
30 // Michael Kelsey 31st January 2019 -- Move f <<
31 // <<
32 // Class Description: <<
33 // <<
34 // Model that knows how to draw the magnetic f 31 // Model that knows how to draw the magnetic field.
35 32
36 #include "G4MagneticFieldModel.hh" 33 #include "G4MagneticFieldModel.hh"
>> 34
>> 35 #include "G4VGraphicsScene.hh"
>> 36 #include "G4TransportationManager.hh"
>> 37 #include "G4FieldManager.hh"
37 #include "G4Field.hh" 38 #include "G4Field.hh"
38 #include "G4Point3D.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
>> 50 G4MagneticFieldModel::~G4MagneticFieldModel () {}
>> 51
>> 52 G4MagneticFieldModel::G4MagneticFieldModel
>> 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
>> 65 void G4MagneticFieldModel::DescribeYourselfTo (G4VGraphicsScene& sceneHandler)
>> 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
>> 91 G4TransportationManager* tMgr =
>> 92 G4TransportationManager::GetTransportationManager();
>> 93 assert(tMgr);
>> 94 G4Navigator* navigator = tMgr->GetNavigatorForTracking();
>> 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.
39 143
>> 144 // Working space for GetFieldValue.
>> 145 double position_time[4] = {0,0,0,0};
>> 146 double result[nResults];
40 147
41 // Return magnetic field vector for display << 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.
42 152
43 void G4MagneticFieldModel:: << 153 // Get field values and ascertain maximum field.
44 GetFieldAtLocation(const G4Field* field, const << 154 G4double BFieldMagnitudeMax = -std::numeric_limits<G4double>::max();
45 G4double time, G4Point3D& result) const << 155 for (G4int i = 0; i < nXSamples; i++) {
46 if (!field) return; // No action if no f << 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 }
47 223
48 G4double xyzt[4] = { position.x(), position. << 224 if (BFieldMagnitudeMax <= 0) {
49 G4double BEvals[6] = {0.}; // Field retur << 225 G4cout
50 field->GetFieldValue(xyzt, BEvals); << 226 << "No field in this scene."
>> 227 << G4endl;
>> 228 return;
>> 229 }
51 230
52 result.set(BEvals[0], BEvals[1], BEvals[2]); << 231 if (fRepresentation == Representation::lightArrow) sceneHandler.BeginPrimitives();
53 return; << 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();
54 } 287 }
55 288