Interpolation/src/TrapezoidalCylindricalMFGrid.cc

0001 #include "TrapezoidalCylindricalMFGrid.h"
0002 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0003 #include "LinearGridInterpolator3D.h"
0004 #include "MagneticField/VolumeGeometry/interface/MagExceptions.h"
0005
0006 #include <iostream>
0007
0008 using namespace std;
0009
0010 TrapezoidalCylindricalMFGrid::TrapezoidalCylindricalMFGrid(binary_ifstream& inFile,
0011                                                            const GloballyPositioned<float>& vol)
0012     : MFGrid3D(vol) {
0013   // The parameters read from the data files are given in global coordinates.
0014   // In version 85l, local frame has the same orientation of global frame for the reference
0015   // volume, i.e. the r.f. transformation is only a translation.
0016   // There is therefore no need to convert the field values to local coordinates.
0017   // Check this assumption:
0018   GlobalVector localXDir(frame().toGlobal(LocalVector(1, 0, 0)));
0019   GlobalVector localYDir(frame().toGlobal(LocalVector(0, 1, 0)));
0020
0021   if (localXDir.dot(GlobalVector(1, 0, 0)) > 0.999999 && localYDir.dot(GlobalVector(0, 1, 0)) > 0.999999) {
0022     // "null" rotation - requires no conversion...
0023   } else {
0024     cout << "ERROR: TrapezoidalCylindricalMFGrid: unexpected orientation: x: " << localXDir << " y: " << localYDir
0025          << endl;
0026   }
0027
0028   int n1, n2, n3;
0029   inFile >> n1 >> n2 >> n3;
0030   double xref, yref, zref;
0031   inFile >> xref >> yref >> zref;
0032   double stepx, stepy, stepz;
0033   inFile >> stepx >> stepy >> stepz;
0034
0035   double BasicDistance1[3][3];  // linear step
0036   double BasicDistance2[3][3];  // linear offset
0037   bool easya, easyb, easyc;
0038
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 >> easya >> easyb >> easyc;
0046
0047   vector<BVector> fieldValues;
0048   float Bx, By, Bz;
0049   int nLines = n1 * n2 * n3;
0050   fieldValues.reserve(nLines);
0051   for (int iLine = 0; iLine < nLines; ++iLine) {
0052     inFile >> Bx >> By >> Bz;
0053     fieldValues.push_back(BVector(Bx, By, Bz));
0054   }
0055   // check completeness
0056   string lastEntry;
0057   inFile >> lastEntry;
0058   if (lastEntry != "complete") {
0059     cout << "ERROR during file reading: file is not complete" << endl;
0060   }
0061
0062 #ifdef DEBUG_GRID
0063   cout << "easya " << easya << " easyb " << easyb << " easyc " << easyc << endl;
0064 #endif
0065
0066   if (!easyb || !easyc) {
0067     throw MagGeometryError("TrapezoidalCartesianMFGrid only implemented for first coordinate");
0068   }
0069
0070 #ifdef DEBUG_GRID
0071   cout << "Grid reference point in grid system: " << xref << "," << yref << "," << zref << endl;
0072   cout << "steps " << stepx << "," << stepy << "," << stepz << endl;
0073   cout << "ns " << n1 << "," << n2 << "," << n3 << endl;
0074
0075   for (int i = 0; i < 3; ++i)
0076     for (int j = 0; j < 3; ++j) {
0077       cout << "BasicDistance1[" << i << "][" << j << "] = " << BasicDistance1[i][j] << "BasicDistance2[" << i << "]["
0078            << j << "] = " << BasicDistance2[i][j] << endl;
0079     }
0080 #endif
0081
0082   // the "not easy" coordinate is x
0083   double a = stepx * (n1 - 1);
0084   double b = a + BasicDistance1[0][1] * (n2 - 1) * (n1 - 1) + BasicDistance1[0][2] * (n3 - 1) * (n1 - 1);
0085   //  double h = stepy * (n2-1);
0086   double h = stepz * (n3 - 1);
0087   double delta = -BasicDistance2[0][1] * (n2 - 1) - BasicDistance2[0][2] * (n3 - 1);
0088
0089 #ifdef DEBUG_GRID
0090   cout << "Trapeze size (a,b,h) = " << a << "," << b << "," << h << endl;
0091 #endif
0092
0093   GlobalPoint grefp(GlobalPoint::Cylindrical(xref, Geom::pi() - yref, zref));
0094   LocalPoint lrefp = frame().toLocal(grefp);
0095
0096 #ifdef DEBUG_GRID
0097   cout << "Global origin " << grefp << endl;
0098   cout << "Local origin  " << lrefp << endl;
0099 #endif
0100
0101   double baMinus1 = BasicDistance1[0][2] * (n3 - 1) / stepx;
0102   if (std::abs(baMinus1) > 0.000001) {
0103     double b_over_a = 1 + baMinus1;
0104     double a1 = std::abs(baMinus1) > 0.000001 ? delta / baMinus1 : a / 2;
0105 #ifdef DEBUG_GRID
0106     cout << "a1 = " << a1 << endl;
0107 #endif
0108
0109     // transform reference point to grid frame
0110     double x0 = lrefp.perp() + a1;
0111     double y0 = lrefp.z() + h / 2.;
0112     mapping_ = Trapezoid2RectangleMappingX(x0, y0, b_over_a, h);
0113   } else {  // parallelogram
0114     mapping_ = Trapezoid2RectangleMappingX(0, 0, delta / h);
0115   }
0116   double xrec, yrec;
0117   mapping_.rectangle(lrefp.perp(), lrefp.z(), xrec, yrec);
0118
0119   Grid1D gridX(xrec, xrec + (a + b) / 2., n1);
0120   Grid1D gridY(yref, yref + stepy * (n2 - 1), n2);
0121   Grid1D gridZ(yrec, yrec + h, n3);
0122   grid_ = GridType(gridX, gridY, gridZ, fieldValues);
0123
0124   // Activate/deactivate timers
0125   //   static SimpleConfigurable<bool> timerOn(false,"MFGrid:timing");
0126   //   (*TimingReport::current()).switchOn("MagneticFieldProvider::valueInTesla(TrapezoidalCylindricalMFGrid)",timerOn);
0127 }
0128
0129 void TrapezoidalCylindricalMFGrid::dump() const {}
0130
0131 MFGrid::LocalVector TrapezoidalCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0132   //   static TimingReport::Item & timer= (*TimingReport::current())["MagneticFieldProvider::valueInTesla(TrapezoidalCylindricalMFGrid)"];
0133   //   TimeMe t(timer,false);
0134
0135   LinearGridInterpolator3D interpol(grid_);
0136   double a, b, c;
0137   toGridFrame(p, a, b, c);
0138   GlobalVector gv(interpol.interpolate(a, b, c));  // grid in global frame
0139   return frame().toLocal(gv);                      // must return a local vector
0140 }
0141
0142 void TrapezoidalCylindricalMFGrid::toGridFrame(const LocalPoint& p, double& a, double& b, double& c) const {
0143   mapping_.rectangle(p.perp(), p.z(), a, c);
0144   // FIXME: "OLD" convention of phi.
0145   //  b = Geom::pi() - p.phi();
0146   b = p.phi();
0147 }
0148
0149 MFGrid::LocalPoint TrapezoidalCylindricalMFGrid::fromGridFrame(double a, double b, double c) const {
0150   double rtrap, ztrap;
0151   mapping_.trapezoid(a, c, rtrap, ztrap);
0152   // FIXME: "OLD" convention of phi.
0153   //  return LocalPoint(LocalPoint::Cylindrical(rtrap, Geom::pi() - b, ztrap));
0154   return LocalPoint(LocalPoint::Cylindrical(rtrap, b, ztrap));
0155 }