Core/include/eASTMagneticField.hh

0001 // ********************************************************************
0002 //
0003 // eASTMagneticField.hh
0004 //   Header file of eAST Magnetic Field class
0005 //
0006 // History
0007 //   June 22nd, 2021 : first implementation
0008 //
0009 // ********************************************************************
0010
0011 #ifndef eASTMagneticField_h
0012 #define eASTMagneticField_h 1
0013
0014 #include "G4MagneticField.hh"
0015 #include "G4GenericMessenger.hh"
0016 #include "G4SystemOfUnits.hh"
0017 #include "G4UnitsTable.hh"
0018
0019 #include <vector>
0020 #include <array>
0021
0022 class G4EquationOfMotion;
0023 class G4MagIntegratorStepper;
0024 class G4ChordFinder;
0025 class G4FieldManager;
0026 class G4PropagatorInField;
0027
0028 class eASTMagneticFieldMap {
0029
0030   public:
0031     eASTMagneticFieldMap(const G4String& name);
0032     ~eASTMagneticFieldMap() = default;
0033
0034   private:
0035
0036     // Static maps to indices of neighboring points
0037     static const int kLinearMap[8][3];
0038     static const int kCubicMap[64][3];
0039
0040   private:
0041
0042     // Field name
0043     G4String fName{""};
0044
0045   public:
0046
0047     // Load field map from filename
0048     bool Load(const G4String& filename);
0049
0050     // Get field value at point
0051     void GetFieldValue(const G4double point[4], G4double *field) const {
0052       field[0] = 0;
0053       field[1] = 0;
0054       field[2] = 0;
0055       AddFieldValue(point, field);
0056     };
0057
0058     // Add field value at point
0059     void AddFieldValue(const G4double point[4], G4double *field) const;
0060
0061     // Print field value at point
0062     void PrintFieldValue(const G4ThreeVector& point);
0063
0064   private:
0065
0066     // Field messenger
0067     G4GenericMessenger fFieldMapMessenger;
0068
0069   private:
0070
0071     // File format
0072     enum EFileFormat {kFileFormatUndefined, kFileFormatRZPhi, kFileFormatRZ};
0073     EFileFormat fFileFormat{kFileFormatUndefined};
0074
0075   private:
0076
0077     // Map units and extent
0078     double fFieldUnit{CLHEP::tesla};
0079     std::array<double,3> fGridUnit{CLHEP::cm, CLHEP::cm, CLHEP::deg};
0080     std::array<std::tuple<double,double,double>,3> fGridExtent; // min, max, spacing
0081     std::array<unsigned int,3> fGridSize{0,0,0};
0082
0083   public:
0084
0085     // Map units and extent getters and setters
0086     void SetFieldUnit(const G4String& unit) {
0087       fFieldUnit = G4UnitDefinition::GetValueOf(unit);
0088     }
0089     void SetGridUnit(unsigned int i, const G4String& unit) {
0090       if (i < fGridExtent.size()) {
0091         fGridUnit[i] = G4UnitDefinition::GetValueOf(unit);
0092       } else {
0093         G4cerr << "ERROR: cannot set magnetic field grid unit" << G4endl;
0094       }
0095     }
0096     void SetGridExtent(unsigned int i, const G4ThreeVector& extent) {
0097       if (i < fGridExtent.size() && extent.x() < extent.y() && extent.z() > 0.0) {
0098         fGridExtent[i] = std::make_tuple(extent.x(), extent.y(), extent.z());
0099       } else {
0100         G4cerr << "ERROR: cannot set magnetic field grid extent" << G4endl;
0101       }
0102     }
0103
0104   private:
0105
0106     // Map data
0107     std::vector<std::vector<std::vector<std::vector<double>>>> fMap;
0108
0109   private:
0110
0111     // Interpolation type to allow for linear and cubic interpolation
0112     enum EInterpolationType {kLinear, kCubic};
0113     EInterpolationType fInterpolationType{kLinear};
0114
0115   public:
0116
0117     // Interpolation type setters and getters
0118     void SetInterpolationType(const G4String& type) {
0119       if (type == "linear") fInterpolationType = kLinear;
0120       if (type == "cubic") fInterpolationType = kCubic;
0121     }
0122     EInterpolationType GetInterpolationType() const {
0123       return fInterpolationType;
0124     }
0125
0126   private:
0127
0128     // Linear interpolation functions
0129     double _linearInterpolate(const double p[2 << 0], double x) const {
0130       return p[0] + x * (p[1] - p[0]);
0131     }
0132     double _bilinearInterpolate(const double p[2 << 1], const double x[2]) const {
0133       double c[2];
0134       c[0] = _linearInterpolate(&(p[0]), x[1]);
0135       c[1] = _linearInterpolate(&(p[2]), x[1]);
0136       return _linearInterpolate(c, x[0]);
0137     }
0138     double _trilinearInterpolate(const double p[2 << 2], const double x[3]) const {
0139       double c[2];
0140       c[0] = _bilinearInterpolate(&(p[0]), &(x[1]));
0141       c[1] = _bilinearInterpolate(&(p[4]), &(x[1]));
0142       return _linearInterpolate(c, x[0]);
0143     }
0144
0145     // Cubic interpolation functions
0146     double _cubicInterpolate(const double p[4 << 0], double x) const {
0147       return p[1] + 0.5 * x * (p[2] - p[0] +
0148                                x * (2. * p[0] - 5. * p[1] + 4. * p[2] - p[3] +
0149                                     x * (3. * (p[1] - p[2]) + p[3] - p[0])));
0150     }
0151     double _bicubicInterpolate(const double p[4 << 1], const double x[2]) const {
0152       double c[4];
0153       c[0] = _cubicInterpolate(&(p[0]),  x[1]);
0154       c[1] = _cubicInterpolate(&(p[4]),  x[1]);
0155       c[2] = _cubicInterpolate(&(p[8]),  x[1]);
0156       c[3] = _cubicInterpolate(&(p[12]), x[1]);
0157       return _cubicInterpolate(c, x[0]);
0158     }
0159     double _tricubicInterpolate(const double p[4 << 2], const double x[3]) const {
0160       double c[4];
0161       c[0] = _bicubicInterpolate(&(p[0]),  &(x[1]));
0162       c[1] = _bicubicInterpolate(&(p[16]), &(x[1]));
0163       c[2] = _bicubicInterpolate(&(p[32]), &(x[1]));
0164       c[3] = _bicubicInterpolate(&(p[48]), &(x[1]));
0165       return _cubicInterpolate(c, x[0]);
0166     }
0167
0168 };
0169
0170 class eASTMagneticField : public G4MagneticField {
0171
0172   public:
0173
0174     eASTMagneticField();
0175     ~eASTMagneticField() = default;
0176
0177     // Activate in ConstuctSDandField
0178     void Activate();
0179
0180     // Create new field
0181     void CreateField(const G4String& name);
0182
0183     // Get field value at point
0184     void GetFieldValue(const G4double point[4], G4double *field) const {
0185       field[0] = 0;
0186       field[1] = 0;
0187       field[2] = 0;
0188       for (auto &map: fMaps) {
0189         map.AddFieldValue(point, field);
0190       }
0191     };
0192
0193     // Print field value at point
0194     void PrintFieldValue(const G4ThreeVector& point) {
0195       for (auto &map: fMaps) {
0196         map.PrintFieldValue(point);
0197       }
0198     };
0199
0200   private:
0201
0202     // Field messenger
0203     G4GenericMessenger fFieldMessenger;
0204
0205   private:
0206
0207     // Field maps
0208     std::vector<eASTMagneticFieldMap> fMaps;
0209
0210   private:
0211
0212     G4double fMinStep{0.01*mm};
0213     G4double fDeltaChord{3.0*mm};
0214     G4double fDeltaOneStep{0.01*mm};
0215     G4double fDeltaIntersection{0.1*mm};
0216     G4double fEpsMin{1.0e-5*mm};
0217     G4double fEpsMax{1.0e-4*mm};
0218
0219     G4EquationOfMotion*     fEquation{nullptr};
0220     G4int                   fEquationDoF{0};
0221     G4FieldManager*         fFieldManager{nullptr};
0222     G4PropagatorInField*    fFieldPropagator{nullptr};
0223     G4MagIntegratorStepper* fStepper{nullptr};
0224     G4ChordFinder*          fChordFinder{nullptr};
0225
0226 };
0227
0228 #endif