Interpolation/src/MagneticFieldGrid.cc

0001 // include header for MagneticFieldGrid (regular + extension for some trapezoids)
0002 #include "MagneticFieldGrid.h"
0003 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0004 #include <cassert>
0005
0006 using namespace std;
0007
0008 void MagneticFieldGrid::load(const string &name) {
0009   magneticfield::interpolation::binary_ifstream inFile(name);
0010   inFile >> GridType;
0011   // reading the header
0012   switch (GridType) {
0013     case 1:
0014       inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
0015       inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
0016       inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
0017       break;
0018     case 2:
0019       inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
0020       inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
0021       inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
0022       inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
0023       inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
0024       inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
0025       inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
0026       inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
0027       inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
0028       inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
0029       break;
0030     case 3:
0031       inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
0032       inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
0033       inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
0034       break;
0035     case 4:
0036       inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
0037       inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
0038       inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
0039       inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
0040       inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
0041       inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
0042       inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
0043       inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
0044       inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
0045       inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
0046       break;
0047     case 5:
0048       inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
0049       inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
0050       inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
0051       inFile >> RParAsFunOfPhi[0] >> RParAsFunOfPhi[1] >> RParAsFunOfPhi[2] >> RParAsFunOfPhi[3];
0052       break;
0053     default:
0054       assert(0);  //this is a bug
0055   }
0056   //reading the field
0057   float Bx, By, Bz;
0058   BVector FieldEntry;
0059   int nLines = NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2];
0060   FieldValues.reserve(nLines);
0061
0062   for (int iLine = 0; iLine < nLines; ++iLine) {
0063     inFile >> Bx >> By >> Bz;
0064     FieldEntry.putB3(Bx, By, Bz);
0065     FieldValues.push_back(FieldEntry);
0066   }
0067   // check completeness and close file
0068   string lastEntry;
0069   inFile >> lastEntry;
0070   inFile.close();
0071   if (lastEntry != "complete") {
0072     GridType = 0;
0073     cout << "error during file reading: file is not complete" << endl;
0074   }
0075   return;
0076 }
0077
0078 int MagneticFieldGrid::gridType() {
0079   int type = GridType;
0080   bool text = false;
0081   if (text) {
0082     if (type == 0)
0083       cout << "  grid type = " << type << "  -->  not determined" << endl;
0084     if (type == 1)
0085       cout << "  grid type = " << type << "  -->  (x,y,z) cube" << endl;
0086     if (type == 2)
0087       cout << "  grid type = " << type << "  -->  (x,y,z) trapezoid" << endl;
0088     if (type == 3)
0089       cout << "  grid type = " << type << "  -->  (r,phi,z) cube" << endl;
0090     if (type == 4)
0091       cout << "  grid type = " << type << "  -->  (r,phi,z) trapezoid" << endl;
0092     if (type == 5)
0093       cout << "  grid type = " << type << "  -->  (r,phi,z) 1/sin(phi)" << endl;
0094   }
0095   return type;
0096 }
0097
0098 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz) {
0099   double dB[3] = {0., 0., 0.};
0100   // define interpolation object
0101   VectorFieldInterpolation MagInterpol;
0102   // calculate indices for "CellPoint000"
0103   int index[3];
0104   putCoordGetInd(X1, X2, X3, index[0], index[1], index[2]);
0105   int index0[3] = {0, 0, 0};
0106   int index1[3] = {0, 0, 0};
0107   for (int i = 0; i < 3; ++i) {
0108     if (NumberOfPoints[i] > 1) {
0109       index0[i] = max(0, index[i]);
0110       if (index0[i] > NumberOfPoints[i] - 2)
0111         index0[i] = NumberOfPoints[i] - 2;
0112       index1[i] = max(1, index[i] + 1);
0113       if (index1[i] > NumberOfPoints[i] - 1)
0114         index1[i] = NumberOfPoints[i] - 1;
0115     }
0116   }
0117   double tmpX[3];
0118   float tmpB[3];
0119   // define the corners of interpolation volume
0120   // FIXME: should not unpack the arrays to then repack them as first thing.
0121   putIndicesGetB(index0[0], index0[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
0122   putIndGetCoord(index0[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
0123   MagInterpol.defineCellPoint000(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0124   putIndicesGetB(index1[0], index0[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
0125   putIndGetCoord(index1[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
0126   MagInterpol.defineCellPoint100(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0127   putIndicesGetB(index0[0], index1[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
0128   putIndGetCoord(index0[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
0129   MagInterpol.defineCellPoint010(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0130   putIndicesGetB(index1[0], index1[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
0131   putIndGetCoord(index1[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
0132   MagInterpol.defineCellPoint110(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0133   putIndicesGetB(index0[0], index0[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
0134   putIndGetCoord(index0[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
0135   MagInterpol.defineCellPoint001(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0136   putIndicesGetB(index1[0], index0[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
0137   putIndGetCoord(index1[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
0138   MagInterpol.defineCellPoint101(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0139   putIndicesGetB(index0[0], index1[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
0140   putIndGetCoord(index0[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
0141   MagInterpol.defineCellPoint011(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0142   putIndicesGetB(index1[0], index1[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
0143   putIndGetCoord(index1[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
0144   MagInterpol.defineCellPoint111(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
0145   // interpolate
0146   MagInterpol.putSCoordGetVField(X1, X2, X3, dB[0], dB[1], dB[2]);
0147   Bx = float(dB[0]);
0148   By = float(dB[1]);
0149   Bz = float(dB[2]);
0150   return;
0151 }
0152
0153 // FIXME: Signature should be:
0154 //
0155 //     void MagneticFieldGrid::putCoordGetInd(double *pos, int *index)
0156 //
0157 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3) {
0158   double pnt[3] = {X1, X2, X3};
0159   int index[3];
0160   switch (GridType) {
0161     case 1: {
0162       for (int i = 0; i < 3; ++i) {
0163         index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
0164       }
0165       break;
0166     }
0167     case 2: {
0168       // FIXME: Should use else!
0169       for (int i = 0; i < 3; ++i) {
0170         if (EasyCoordinate[i]) {
0171           index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
0172         } else
0173           index[i] = 0;  //computed below
0174       }
0175       for (int i = 0; i < 3; ++i) {
0176         if (!EasyCoordinate[i]) {
0177           double stepSize = BasicDistance0[i];
0178           double offset = 0.0;
0179           for (int j = 0; j < 3; ++j) {
0180             stepSize += BasicDistance1[i][j] * index[j];
0181             offset += BasicDistance2[i][j] * index[j];
0182           }
0183           index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
0184         }
0185       }
0186       break;
0187     }
0188     case 3: {
0189       for (int i = 0; i < 3; ++i) {
0190         index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
0191       }
0192       break;
0193     }
0194     case 4: {
0195       // FIXME: should use else!
0196       for (int i = 0; i < 3; ++i) {
0197         if (EasyCoordinate[i]) {
0198           index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
0199         } else
0200           index[i] = 0;  //computed below
0201       }
0202       for (int i = 0; i < 3; ++i) {
0203         if (!EasyCoordinate[i]) {
0204           double stepSize = BasicDistance0[i];
0205           double offset = 0.0;
0206           for (int j = 0; j < 3; ++j) {
0207             stepSize += BasicDistance1[i][j] * index[j];
0208             offset += BasicDistance2[i][j] * index[j];
0209           }
0210           index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
0211         }
0212       }
0213       break;
0214     }
0215     case 5: {
0216       double sinPhi = sin(pnt[1]);
0217       double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
0218       stepSize = stepSize / (NumberOfPoints[0] - 1);
0219       double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
0220       index[0] = int((pnt[0] - startingPoint) / stepSize);
0221       index[1] = int((pnt[1] - ReferencePoint[1]) / BasicDistance0[1]);
0222       index[2] = int((pnt[2] - ReferencePoint[2]) / BasicDistance0[2]);
0223       break;
0224     }
0225     default:
0226       assert(0);  //shouldn't be here
0227   }
0228   Index1 = index[0];
0229   Index2 = index[1];
0230   Index3 = index[2];
0231   return;
0232 }
0233
0234 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz) {
0235   BVector FieldEntry;
0236   FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
0237   Bx = FieldEntry.bx();
0238   By = FieldEntry.by();
0239   Bz = FieldEntry.bz();
0240   return;
0241 }
0242
0243 // FIXME: same as above.
0244 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3) {
0245   int index[3] = {Index1, Index2, Index3};
0246   double pnt[3];
0247   switch (GridType) {
0248     case 1: {
0249       for (int i = 0; i < 3; ++i) {
0250         pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
0251       }
0252       break;
0253     }
0254     case 2: {
0255       for (int i = 0; i < 3; ++i) {
0256         if (EasyCoordinate[i]) {
0257           pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
0258         } else {
0259           double stepSize = BasicDistance0[i];
0260           double offset = 0.0;
0261           for (int j = 0; j < 3; ++j) {
0262             stepSize += BasicDistance1[i][j] * index[j];
0263             offset += BasicDistance2[i][j] * index[j];
0264           }
0265           pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
0266         }
0267       }
0268       break;
0269     }
0270     case 3: {
0271       for (int i = 0; i < 3; ++i) {
0272         pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
0273       }
0274       break;
0275     }
0276     case 4: {
0277       for (int i = 0; i < 3; ++i) {
0278         if (EasyCoordinate[i]) {
0279           pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
0280         } else {
0281           double stepSize = BasicDistance0[i];
0282           double offset = 0.0;
0283           for (int j = 0; j < 3; ++j) {
0284             stepSize += BasicDistance1[i][j] * index[j];
0285             offset += BasicDistance2[i][j] * index[j];
0286           }
0287           pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
0288         }
0289       }
0290       break;
0291     }
0292     case 5: {
0293       pnt[2] = ReferencePoint[2] + BasicDistance0[2] * index[2];
0294       pnt[1] = ReferencePoint[1] + BasicDistance0[1] * index[1];
0295       double sinPhi = sin(pnt[1]);
0296       double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
0297       stepSize = stepSize / (NumberOfPoints[0] - 1);
0298       double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
0299       pnt[0] = startingPoint + stepSize * index[0];
0300       break;
0301     }
0302     default:
0303       assert(0);  //bug if make it here
0304   }
0305   X1 = pnt[0];
0306   X2 = pnt[1];
0307   X3 = pnt[2];
0308   return;
0309 }
0310
0311 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3) {
0312   return Index1 * NumberOfPoints[1] * NumberOfPoints[2] + Index2 * NumberOfPoints[2] + Index3;
0313 }
0314
0315 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz) {
0316   B3[0] = Bx;
0317   B3[1] = By;
0318   B3[2] = Bz;
0319   return;
0320 }
0321
0322 float MagneticFieldGrid::BVector::bx() { return B3[0]; }
0323
0324 float MagneticFieldGrid::BVector::by() { return B3[1]; }
0325
0326 float MagneticFieldGrid::BVector::bz() { return B3[2]; }