test/BinaryTablesGeneration/prepareMagneticFieldGrid.h

0001 #ifndef prepareMagneticFieldGrid_H
0002 #define prepareMagneticFieldGrid_H
0003
0004 /** \class prepareMagneticFieldGrid
0005  *
0006  * read coordinates and magnetic field values from an ASCII file
0007  *
0008  * The file given as inmust be in the form:
0009  *   *-xyz-* (Cartesian)
0010  *   *-rpz-* (Cylindrical)
0011  * -remark: the units of the ASCII file are unknown to this class
0012  *
0013  * Inputs are assumed to be in local coordinates if rotateFromSector=0.
0014  * Otherwise, inputs are assumed to be in global coordinates, and are converted
0015  * to the local reference frame of the sector specified in rotateFromSector.
0016  *
0017  * determine the structure of the coordinate points (e.g. trapezoid)
0018  * and store the information on a grid like structure for fast access
0019  * -remark: one variable may depend linearly on the others
0020  * -remark: points do not need to be ordered in ASCII file
0021  *
0022  * additional functions either translate indices <-> coordinates,
0023  * transfer data, or activate the interpolation between grid points
0024  *
0025  * \author : <Volker.Drollinger@cern.ch>, updated N. Amapane 04/2008, 2102, 2013
0026  *
0027  *
0028  */
0029
0030 // used libs
0031 #include <vector>
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <string>
0035 #include <iostream>
0036 #include <fstream>
0037
0038 #define PRINT true
0039 #define EPSILON 1e-4
0040
0041 class prepareMagneticFieldGrid {
0042 public:
0043   // constructor
0044   prepareMagneticFieldGrid(int rotateFromSector = 0) : rotateSector(rotateFromSector) {
0045     GridType = 0;
0046     for (int i = 0; i < 3; ++i) {
0047       NumberOfPoints[i] = 0;
0048     };
0049     for (int i = 0; i < 3; ++i) {
0050       ReferencePoint[i] = 0.;
0051     };
0052     for (int i = 0; i < 3; ++i) {
0053       BasicDistance0[i] = 0.;
0054     };
0055     for (int i = 0; i < 3; ++i) {
0056       for (int j = 0; j < 3; ++j) {
0057         BasicDistance1[i][j] = 0.;
0058       };
0059     };
0060     for (int i = 0; i < 3; ++i) {
0061       for (int j = 0; j < 3; ++j) {
0062         BasicDistance2[i][j] = 0.;
0063       };
0064     };
0065     for (int i = 0; i < 4; ++i) {
0066       RParAsFunOfPhi[i] = 0.;
0067     };
0068     for (int i = 0; i < 3; ++i) {
0069       EasyCoordinate[i] = false;
0070     };
0071     KnownStructure = false;
0072     XyzCoordinates = false;
0073     RpzCoordinates = false;
0074     sector = unknown;
0075   }
0076   // destructor
0077   ~prepareMagneticFieldGrid() {}
0078
0079 private:
0080   enum masterSector { unknown = 0, one = 1, four = 4 };
0081
0082   class IndexedDoubleVector {
0083   public:
0084     // constructor
0085     IndexedDoubleVector() {}
0086     // destructor
0087     ~IndexedDoubleVector() {}
0088     // less operator (defined for indices: 1st index = highest priority, 2nd index = 2nd highest priority, ...)
0089     bool operator<(const IndexedDoubleVector &) const;
0090
0091   private:
0092     // six component vector in double precision
0093     int I3[3];
0094     double V6[6];
0095
0096   public:
0097     // Accessors
0098     void putI3(int index1, int index2, int index3);
0099     void getI3(int &index1, int &index2, int &index3);
0100     void putV6(double X1, double X2, double X3, double Bx, double By, double Bz);
0101     void getV6(double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz);
0102   };
0103   class SixDPoint {
0104   public:
0105     // constructor
0106     SixDPoint() {}
0107     // destructor
0108     ~SixDPoint() {}
0109
0110   private:
0111     // six component vector in double precision
0112     double P6[6];
0113
0114   public:
0115     // Accessors
0116     void putP6(double X1, double X2, double X3, double Bx, double By, double Bz);
0117     double x1();
0118     double x2();
0119     double x3();
0120     double bx();
0121     double by();
0122     double bz();
0123   };
0124   // definition of grid
0125   // structure
0126   int GridType;
0127   int NumberOfPoints[3];
0128   double ReferencePoint[3];
0129   double BasicDistance0[3];     // constant step
0130   double BasicDistance1[3][3];  // linear step
0131   double BasicDistance2[3][3];  // linear offset
0132   double RParAsFunOfPhi[4];     // R = f(phi) or const. (0,2: const. par. ; 1,3: const./sin(phi))
0133   bool EasyCoordinate[3];
0134   bool KnownStructure;
0135   bool XyzCoordinates;
0136   bool RpzCoordinates;
0137   int rotateSector;
0138   masterSector sector;
0139
0140   // all points (X1,X2,X3,Bx,By,Bz) of one volume
0141   std::vector<SixDPoint> GridData;
0142
0143   // convert original units [m] to new units [cm]
0144   void convertUnits();
0145
0146 public:
0147   /// check, if number of lines corresponds to number of points in space (double counting test)
0148   void countTrueNumberOfPoints(const std::string &name) const;
0149   /// reads the corresponding ASCII file, detects the logic of the points, and saves them on a grid
0150   void fillFromFile(const std::string &name);
0151   /// sames as fillFromFile, but for special cases which are not covered by the standard algorithm
0152   void fillFromFileSpecial(const std::string &name);
0153   /// returns value of GridType (and eventually prints the type + short description)
0154   int gridType();
0155   /// possibility to check existing MagneticFieldGrid point by point (all points)
0156   void validateAllPoints();
0157   /// sames as fillFromFile, but for special cases which are not covered by the standard algorithm
0158   void saveGridToFile(const std::string &outName);
0159
0160   /// indicates, that MagneticFieldGrid is fully operational (for interpolation)
0161   bool isReady();
0162
0163   /// interpolates the magnetic field at input coordinate point and returns field values
0164   void interpolateAtPoint(double X1, double X2, double X3, double &Bx, double &By, double &Bz);
0165
0166   // calculates indices from coordinates
0167   void putCoordGetIndices(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3);
0168   // takes indices and returns coordinates and magnetic field values
0169   void putIndicesGetXAndB(
0170       int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz);
0171   // takes indices, calculates coordinates, and returns coordinates + magnetic fiels values (for testing purposes only)
0172   void putIndCalcXReturnB(
0173       int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz);
0174   // converts three indices into one number (for the vector GridData)
0175   int lineNumber(int Index1, int Index2, int Index3);
0176 };
0177
0178 #endif