Interpolation/src/TrapezoidalCartesianMFGrid.cc

0001 /** \file
0002  *
0003  *  \author T. Todorov - updated N. Amapane (2008)
0004  */
0005
0006 #include "TrapezoidalCartesianMFGrid.h"
0007 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0008 #include "LinearGridInterpolator3D.h"
0009 #include "MagneticField/VolumeGeometry/interface/MagExceptions.h"
0010 #include <iostream>
0011
0012 //#define DEBUG_GRID
0013
0014 using namespace std;
0015
0016 TrapezoidalCartesianMFGrid::TrapezoidalCartesianMFGrid(binary_ifstream& inFile, const GloballyPositioned<float>& vol)
0017     : MFGrid3D(vol), increasingAlongX(false), convertToLocal(true) {
0018   // The parameters read from the data files are given in global coordinates.
0019   // In version 85l, local frame has the same orientation of global frame for the reference
0020   // volume, i.e. the r.f. transformation is only a translation.
0021   // There is therefore no need to convert the field values to local coordinates.
0022
0023   // Check orientation of local reference frame:
0024   GlobalVector localXDir(frame().toGlobal(LocalVector(1, 0, 0)));
0025   GlobalVector localYDir(frame().toGlobal(LocalVector(0, 1, 0)));
0026
0027   if (localXDir.dot(GlobalVector(1, 0, 0)) > 0.999999 && localYDir.dot(GlobalVector(0, 1, 0)) > 0.999999) {
0028     // "null" rotation - requires no conversion...
0029     convertToLocal = false;
0030   } else if (localXDir.dot(GlobalVector(0, 1, 0)) > 0.999999 && localYDir.dot(GlobalVector(1, 0, 0)) > 0.999999) {
0031     // Typical orientation if master volume is in sector 1
0032     convertToLocal = true;
0033   } else {
0034     convertToLocal = true;
0035     // Nothing wrong in principle, but this is not expected
0036     cout << "WARNING: TrapezoidalCartesianMFGrid: unexpected orientation: x: " << localXDir << " y: " << localYDir
0037          << endl;
0038   }
0039
0040   int n1, n2, n3;
0041   inFile >> n1 >> n2 >> n3;
0042   double xref, yref, zref;
0043   inFile >> xref >> yref >> zref;
0044   double step1, step2, step3;
0045   inFile >> step1 >> step2 >> step3;
0046
0047   double BasicDistance1[3][3];  // linear step
0048   double BasicDistance2[3][3];  // linear offset
0049   bool easya, easyb, easyc;
0050
0051   inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
0052   inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
0053   inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
0054   inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
0055   inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
0056   inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
0057   inFile >> easya >> easyb >> easyc;
0058
0059   vector<BVector> fieldValues;
0060   float Bx, By, Bz;
0061   int nLines = n1 * n2 * n3;
0062   fieldValues.reserve(nLines);
0063   for (int iLine = 0; iLine < nLines; ++iLine) {
0064     inFile >> Bx >> By >> Bz;
0065     if (convertToLocal) {
0066       // Preserve double precision!
0067       Vector3DBase<double, LocalTag> lB = frame().toLocal(Vector3DBase<double, GlobalTag>(Bx, By, Bz));
0068       fieldValues.push_back(BVector(lB.x(), lB.y(), lB.z()));
0069
0070     } else {
0071       fieldValues.push_back(BVector(Bx, By, Bz));
0072     }
0073   }
0074   // check completeness
0075   string lastEntry;
0076   inFile >> lastEntry;
0077   if (lastEntry != "complete") {
0078     cout << "ERROR during file reading: file is not complete" << endl;
0079   }
0080
0081   // In version 1103l the reference sector is at phi=0 so that local y is along global X.
0082   // The increasing grid steps for the same volume are then along Y instead than along X.
0083   // To use Trapezoid2RectangleMappingX to map the trapezoidal geometry to a rectangular
0084   // cartesian geometry, we have to exchange (global) X and Y appropriately when constructing
0085   // it.
0086   int nx, ny;
0087   double stepx, stepy, stepz;
0088   double dstep, offset;
0089   if (!easya && easyb && easyc) {
0090     // Increasing grid spacing is on x
0091     increasingAlongX = true;
0092     nx = n1;
0093     ny = n2;
0094     stepx = step1;
0095     stepy = step2;
0096     stepz = step3;
0097     dstep = BasicDistance1[0][1];
0098     offset = BasicDistance2[0][1];
0099   } else if (easya && !easyb && easyc) {
0100     // Increasing grid spacing is on y
0101     increasingAlongX = false;
0102     nx = n2;
0103     ny = n1;
0104     stepx = step2;
0105     stepy = step1;
0106     stepz = -step3;
0107     dstep = BasicDistance1[1][0];
0108     offset = BasicDistance2[1][0];
0109
0110   } else {
0111     // Increasing spacing on z or on > 1 coordinate not supported
0112     throw MagGeometryError("TrapezoidalCartesianMFGrid only implemented for first or second coordinate");
0113   }
0114
0115   double a = stepx * (nx - 1);                 // first base
0116   double b = a + dstep * (ny - 1) * (nx - 1);  // second base
0117   double h = stepy * (ny - 1);                 // height
0118   double delta = -offset * (ny - 1);           // offset between two bases
0119   double baMinus1 = dstep * (ny - 1) / stepx;  // (b*a) - 1
0120
0121   GlobalPoint grefp(xref, yref, zref);
0122   LocalPoint lrefp = frame().toLocal(grefp);
0123
0124   if (fabs(baMinus1) > 0.000001) {
0125     double b_over_a = 1 + baMinus1;
0126     double a1 = delta / baMinus1;
0127
0128 #ifdef DEBUG_GRID
0129     cout << "Trapeze size (a,b,h) = " << a << "," << b << "," << h << endl;
0130     cout << "Global origin " << grefp << endl;
0131     cout << "Local origin  " << lrefp << endl;
0132     cout << "a1 = " << a1 << endl;
0133 #endif
0134
0135     // FIXME ASSUMPTION: here we assume that the local reference frame is oriented with X along
0136     // the direction of where the grid is not uniform. This is the case for all current geometries
0137     double x0 = lrefp.x() + a1;
0138     double y0 = lrefp.y() + h / 2.;
0139     mapping_ = Trapezoid2RectangleMappingX(x0, y0, b_over_a, h);
0140   } else {  // parallelogram
0141     mapping_ = Trapezoid2RectangleMappingX(0, 0, delta / h);
0142   }
0143
0144   // transform reference point to grid frame
0145   double xrec, yrec;
0146   mapping_.rectangle(lrefp.x(), lrefp.y(), xrec, yrec);
0147
0148   Grid1D gridX(xrec, xrec + (a + b) / 2., nx);
0149   Grid1D gridY(yrec, yrec + h, ny);
0150   Grid1D gridZ(lrefp.z(), lrefp.z() + stepz * (n3 - 1), n3);
0151
0152 #ifdef DEBUG_GRID
0153   cout << " GRID X range: local " << gridX.lower() << " - " << gridX.upper()
0154        << " global: " << (frame().toGlobal(LocalPoint(gridX.lower(), 0, 0))).y() << " - "
0155        << (frame().toGlobal(LocalPoint(gridX.upper(), 0, 0))).y() << endl;
0156
0157   cout << " GRID Y range: local " << gridY.lower() << " - " << gridY.upper()
0158        << " global: " << (frame().toGlobal(LocalPoint(0, gridY.lower(), 0))).x() << " - "
0159        << (frame().toGlobal(LocalPoint(0, gridY.upper(), 0))).x() << endl;
0160
0161   cout << " GRID Z range: local " << gridZ.lower() << " - " << gridZ.upper()
0162        << " global: " << (frame().toGlobal(LocalPoint(0, 0, gridZ.lower()))).z() << " "
0163        << (frame().toGlobal(LocalPoint(0, 0, gridZ.upper()))).z() << endl;
0164 #endif
0165
0166   if (increasingAlongX) {
0167     grid_ = GridType(gridX, gridY, gridZ, fieldValues);
0168   } else {
0169     // The reason why gridY and gridX have to be exchanged is because Grid3D::index(i,j,k)
0170     // assumes a specific order for the fieldValues, and we cannot rearrange this vector.
0171     // Given that we exchange grids, we will have to exchange the outpouts of mapping_rectangle()
0172     // and the inputs of mapping_.trapezoid() in the following...
0173     grid_ = GridType(gridY, gridX, gridZ, fieldValues);
0174   }
0175
0176 #ifdef DEBUG_GRID
0177   dump();
0178 #endif
0179 }
0180
0181 void TrapezoidalCartesianMFGrid::dump() const {
0182   cout << endl << "Dump of TrapezoidalCartesianMFGrid" << endl;
0183   //   cout << "Number of points from file   "
0184   //        << n1 << " " << n2 << " " << n3 << endl;
0185   cout << "Number of points from Grid1D " << grid_.grida().nodes() << " " << grid_.gridb().nodes() << " "
0186        << grid_.gridc().nodes() << endl;
0187
0188   //   cout << "Reference Point from file   "
0189   //        << xref << " " << yref << " " << zref << endl;
0190   cout << "Reference Point from Grid1D " << grid_.grida().lower() << " " << grid_.gridb().lower() << " "
0191        << grid_.gridc().lower() << endl;
0192
0193   //   cout << "Basic Distance from file   "
0194   //        <<  stepx << " " << stepy << " " << stepz << endl;
0195   cout << "Basic Distance from Grid1D " << grid_.grida().step() << " " << grid_.gridb().step() << " "
0196        << grid_.gridc().step() << endl;
0197
0198   cout << "Dumping " << grid_.data().size() << " field values " << endl;
0199   // grid_.dump();
0200
0201   // Dump ALL grid points and values
0202   // CAVEAT: if convertToLocal = true in the ctor, points have been converted to LOCAL
0203   // coordinates. To match those from .table files they have to be converted back to global
0204   //   for (int j=0; j < grid_.gridb().nodes(); j++) {
0205   //     for (int i=0; i < grid_.grida().nodes(); i++) {
0206   //       for (int k=0; k < grid_.gridc().nodes(); k++) {
0207   //    cout << i << " " << j << " " << k << " "
0208   //         << frame().toGlobal(LocalPoint(nodePosition(i,j,k))) << " "
0209   //         << nodeValue(i,j,k) << endl;
0210   //       }
0211   //     }
0212   //   }
0213 }
0214
0215 MFGrid::LocalVector TrapezoidalCartesianMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0216   double xrec, yrec;
0217   mapping_.rectangle(p.x(), p.y(), xrec, yrec);
0218
0219   // cout << "TrapezoidalCartesianMFGrid::valueInTesla at local point " << p << endl;
0220   //   cout << p.x() << " " << p.y()
0221   //        << " transformed to grid frame: " << xrec << " " << yrec << endl;
0222
0223   LinearGridInterpolator3D interpol(grid_);
0224
0225   if (!increasingAlongX) {
0226     std::swap(xrec, yrec);
0227     // B values are already converted to local coord!!! otherwise we should:
0228     // GridType::ValueType value = interpol.interpolate( xrec, yrec, p.z());
0229     // return LocalVector(value.y(),value.x(),-value.z());
0230   }
0231
0232   return LocalVector(interpol.interpolate(xrec, yrec, p.z()));
0233 }
0234
0235 void TrapezoidalCartesianMFGrid::toGridFrame(const LocalPoint& p, double& a, double& b, double& c) const {
0236   if (increasingAlongX) {
0237     mapping_.rectangle(p.x(), p.y(), a, b);
0238   } else {
0239     mapping_.rectangle(p.x(), p.y(), b, a);
0240   }
0241
0242   c = p.z();
0243 }
0244
0245 MFGrid::LocalPoint TrapezoidalCartesianMFGrid::fromGridFrame(double a, double b, double c) const {
0246   double xtrap, ytrap;
0247   if (increasingAlongX) {
0248     mapping_.trapezoid(a, b, xtrap, ytrap);
0249   } else {
0250     mapping_.trapezoid(b, a, xtrap, ytrap);
0251   }
0252   return LocalPoint(xtrap, ytrap, c);
0253 }