Geant4 Cross Reference |
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 // File: CCalMagneticField.cc 28 // Description: User Field class implementation. 29 /////////////////////////////////////////////////////////////////////////////// 30 #include <fstream> 31 32 #include "CCalMagneticField.hh" 33 #include "CCalutils.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4FieldManager.hh" 36 37 //#define ddebug 38 //#define debug 39 40 //Constructor and destructor: 41 42 CCalMagneticField::CCalMagneticField(const G4String &filename) : 43 fval(0), pos(0), slope(0), intercept(0) 44 { 45 #ifdef debug 46 fVerbosity = 1; 47 #else 48 fVerbosity = 0; 49 #endif 50 51 //Let's open the file 52 G4cout << " ==> Opening file " << filename << " to read magnetic field..." 53 << G4endl; 54 G4String pathName = std::getenv("CCAL_GLOBALPATH"); 55 std::ifstream is; 56 G4bool ok = openGeomFile(is, pathName, filename); 57 58 if (ok) { 59 findDO(is, G4String("FLDM")); 60 is >> fval >> npts >> xoff; 61 62 if (fVerbosity) 63 G4cout << "Field value " << fval << " # points " << npts << 64 " offset in x " 65 << xoff*mm << G4endl; 66 67 if (npts > 0) { 68 pos = new G4double[npts]; 69 slope = new G4double[npts]; 70 intercept = new G4double[npts]; 71 72 for (G4int i = 0; i < npts; i++) { 73 is >> pos[i] >> slope[i] >> intercept[i]; 74 if (fVerbosity) 75 G4cout << tab << "Position " << i << " " << pos[i] << " Slope " 76 << slope[i] << " Intercept " << intercept[i] << G4endl; 77 } 78 } 79 80 /////////////////////////////////////////////////////////////// 81 // Close the file 82 G4cout << " ==> Closing file " << filename << G4endl; 83 is.close(); 84 } 85 } 86 87 88 CCalMagneticField::~CCalMagneticField() { 89 if (pos) 90 delete[] pos; 91 if (slope) 92 delete[] slope; 93 if (intercept) 94 delete[] intercept; 95 } 96 97 98 // Member functions 99 100 void CCalMagneticField::MagneticField(const G4double x[3], G4double B[3]) const 101 { 102 G4int i=0; 103 for (i=0; i<2; i++) { 104 B[i] = 0*kilogauss; 105 } 106 107 G4double m1=0; 108 G4double c1=0; 109 G4double xnew = x[0]/mm + xoff; 110 if (npts > 0) { 111 for (i=0; i<npts; i++) { 112 if (xnew > pos[i]*mm) { 113 m1 = slope[i]; 114 c1 = intercept[i]; 115 } 116 } 117 } 118 G4double scor = c1 + m*xnew; 119 if (scor < 0.) scor = 0.; 120 if (scor > 1.) scor = 1.0; 121 122 B[2] = scor*fval*kilogauss; 123 if (fVerbosity) 124 { 125 126 G4cout << "Field at x: " << x[0]/mm << "mm (" << xnew << ") = " << 127 B[2]/tesla 128 << "T (m = " << m1 << ", c = " << 129 c1 << ", scale = " << scor << ")" 130 << G4endl; 131 } 132 } 133 134 135 CLHEP::Hep3Vector CCalMagneticField:: 136 MagneticField(const CLHEP::Hep3Vector point) const { 137 138 G4double x[3],B[3]; 139 CLHEP::Hep3Vector v; 140 141 x[0] = point.x(); 142 x[1] = point.y(); 143 x[2] = point.z(); 144 CCalMagneticField::MagneticField(x, B); 145 v.setX(B[0]); 146 v.setY(B[1]); 147 v.setZ(B[2]); 148 return v; 149 } 150 151 152 void CCalMagneticField::GetFieldValue(const G4double x[3], G4double* B) const { 153 CCalMagneticField::MagneticField(x, B); 154 } 155 156