Interpolation/src/RectangularCylindricalMFGrid.cc

0001 #include "RectangularCylindricalMFGrid.h"
0002 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0003 #include "LinearGridInterpolator3D.h"
0004 #include <iostream>
0005
0006 using namespace std;
0007
0008 RectangularCylindricalMFGrid::RectangularCylindricalMFGrid(binary_ifstream& inFile,
0009                                                            const GloballyPositioned<float>& vol)
0010     : MFGrid3D(vol) {
0011   // The parameters read from the data files are given in global coordinates.
0012   // In version 85l, local frame has the same orientation of global frame for the reference
0013   // volume, i.e. the r.f. transformation is only a translation.
0014   // There is therefore no need to convert the field values to local coordinates.
0015   // Check this assumption:
0016   GlobalVector localXDir(frame().toGlobal(LocalVector(1, 0, 0)));
0017   GlobalVector localYDir(frame().toGlobal(LocalVector(0, 1, 0)));
0018
0019   if (localXDir.dot(GlobalVector(1, 0, 0)) > 0.999999 && localYDir.dot(GlobalVector(0, 1, 0)) > 0.999999) {
0020     // "null" rotation - requires no conversion...
0021   } else {
0022     cout << "ERROR: RectangularCylindricalMFGrid: unexpected orientation: x: " << localXDir << " y: " << localYDir
0023          << endl;
0024   }
0025
0026   int n1, n2, n3;
0027   inFile >> n1 >> n2 >> n3;
0028   double xref, yref, zref;
0029   inFile >> xref >> yref >> zref;
0030   double stepx, stepy, stepz;
0031   inFile >> stepx >> stepy >> stepz;
0032
0033   vector<BVector> fieldValues;
0034   float Bx, By, Bz;
0035   int nLines = n1 * n2 * n3;
0036   fieldValues.reserve(nLines);
0037   for (int iLine = 0; iLine < nLines; ++iLine) {
0038     inFile >> Bx >> By >> Bz;
0039     fieldValues.push_back(BVector(Bx, By, Bz));
0040   }
0041   // check completeness
0042   string lastEntry;
0043   inFile >> lastEntry;
0044   if (lastEntry != "complete") {
0045     cout << "ERROR during file reading: file is not complete" << endl;
0046   }
0047
0048   GlobalPoint grefp(GlobalPoint::Cylindrical(xref, yref, zref));
0049   LocalPoint lrefp = frame().toLocal(grefp);
0050
0051 #ifdef DEBUG_GRID
0052   cout << "Grid reference point in grid system: " << xref << "," << yref << "," << zref << endl;
0053   cout << "Grid reference point in global x,y,z: " << grefp << endl;
0054   cout << "Grid reference point in local x,y,z: " << lrefp << endl;
0055   cout << "steps " << stepx << "," << stepy << "," << stepz << endl;
0056 #endif
0057
0058   Grid1D gridX(lrefp.perp(), lrefp.perp() + stepx * (n1 - 1), n1);
0059   //Grid1D gridY( lrefp.phi(), lrefp.phi() + stepy*(n2-1), n2); // wrong: gives zero
0060   Grid1D gridY(yref, yref + stepy * (n2 - 1), n2);
0061   Grid1D gridZ(lrefp.z(), lrefp.z() + stepz * (n3 - 1), n3);
0062
0063   grid_ = GridType(gridX, gridY, gridZ, fieldValues);
0064 }
0065
0066 void RectangularCylindricalMFGrid::dump() const {
0067   cout << endl << "Dump of RectangularCylindricalMFGrid" << endl;
0068   cout << "Number of points from Grid1D " << grid_.grida().nodes() << " " << grid_.gridb().nodes() << " "
0069        << grid_.gridc().nodes() << endl;
0070
0071   cout << "Reference Point from Grid1D " << grid_.grida().lower() << " " << grid_.gridb().lower() << " "
0072        << grid_.gridc().lower() << endl;
0073
0074   cout << "Basic Distance from Grid1D " << grid_.grida().step() << " " << grid_.gridb().step() << " "
0075        << grid_.gridc().step() << endl;
0076
0077   cout << "Dumping " << grid_.data().size() << " field values " << endl;
0078   // grid_.dump();
0079 }
0080
0081 MFGrid::LocalVector RectangularCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0082   const float minimalSignificantR = 1e-6;  // [cm], points below this radius are treated as zero radius
0083   float R = p.perp();
0084   if (R < minimalSignificantR) {
0085     if (grid_.grida().lower() < minimalSignificantR) {
0086       int k = grid_.gridc().index(p.z());
0087       double u = (p.z() - grid_.gridc().node(k)) / grid_.gridc().step();
0088       LocalVector result((1 - u) * grid_(0, 0, k) + u * grid_(0, 0, k + 1));
0089       return result;
0090     }
0091   }
0092
0093   LinearGridInterpolator3D interpol(grid_);
0094   // FIXME: "OLD" convention of phi.
0095   // GridType::ValueType value = interpol( R, Geom::pi() - p.phi(), p.z());
0096   GridType::ReturnType value = interpol.interpolate(R, p.phi(), p.z());
0097   return LocalVector(value);
0098 }
0099
0100 void RectangularCylindricalMFGrid::toGridFrame(const LocalPoint& p, double& a, double& b, double& c) const {
0101   a = p.perp();
0102   // FIXME: "OLD" convention of phi.
0103   //  b = Geom::pi() - p.phi();
0104   b = p.phi();
0105   c = p.z();
0106 }
0107
0108 MFGrid::LocalPoint RectangularCylindricalMFGrid::fromGridFrame(double a, double b, double c) const {
0109   // FIXME: "OLD" convention of phi.
0110   //  return LocalPoint( LocalPoint::Cylindrical(a, Geom::pi() - b, c));
0111   return LocalPoint(LocalPoint::Cylindrical(a, b, c));
0112 }