File indexing completed on 2024-04-06 12:22:32
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
0012
0013
0014
0015
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
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
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
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
0079 }
0080
0081 MFGrid::LocalVector RectangularCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0082 const float minimalSignificantR = 1e-6;
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
0095
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
0103
0104 b = p.phi();
0105 c = p.z();
0106 }
0107
0108 MFGrid::LocalPoint RectangularCylindricalMFGrid::fromGridFrame(double a, double b, double c) const {
0109
0110
0111 return LocalPoint(LocalPoint::Cylindrical(a, b, c));
0112 }