File indexing completed on 2024-04-06 12:22:32
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
0014
0015
0016
0017
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
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];
0036 double BasicDistance2[3][3];
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
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
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
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
0110 double x0 = lrefp.perp() + a1;
0111 double y0 = lrefp.z() + h / 2.;
0112 mapping_ = Trapezoid2RectangleMappingX(x0, y0, b_over_a, h);
0113 } else {
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
0125
0126
0127 }
0128
0129 void TrapezoidalCylindricalMFGrid::dump() const {}
0130
0131 MFGrid::LocalVector TrapezoidalCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0132
0133
0134
0135 LinearGridInterpolator3D interpol(grid_);
0136 double a, b, c;
0137 toGridFrame(p, a, b, c);
0138 GlobalVector gv(interpol.interpolate(a, b, c));
0139 return frame().toLocal(gv);
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
0145
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
0153
0154 return LocalPoint(LocalPoint::Cylindrical(rtrap, b, ztrap));
0155 }