File indexing completed on 2024-04-06 12:22:32
0001 #include "SpecialCylindricalMFGrid.h"
0002 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0003 #include "LinearGridInterpolator3D.h"
0004
0005 #include <iostream>
0006
0007 using namespace std;
0008
0009 SpecialCylindricalMFGrid::SpecialCylindricalMFGrid(binary_ifstream& inFile,
0010 const GloballyPositioned<float>& vol,
0011 int gridType)
0012 : MFGrid3D(vol) {
0013 if (gridType == 5) {
0014 sector1 = false;
0015 } else if (gridType == 6) {
0016 sector1 = true;
0017 } else {
0018 cout << "ERROR wrong SpecialCylindricalMFGrid type " << gridType << endl;
0019 sector1 = false;
0020 }
0021
0022 int n1, n2, n3;
0023 inFile >> n1 >> n2 >> n3;
0024 #ifdef DEBUG_GRID
0025 cout << "n1 " << n1 << " n2 " << n2 << " n3 " << n3 << endl;
0026 #endif
0027 double xref, yref, zref;
0028 inFile >> xref >> yref >> zref;
0029 double stepx, stepy, stepz;
0030 inFile >> stepx >> stepy >> stepz;
0031
0032 double RParAsFunOfPhi[4];
0033 inFile >> RParAsFunOfPhi[0] >> RParAsFunOfPhi[1] >> RParAsFunOfPhi[2] >> RParAsFunOfPhi[3];
0034
0035 vector<BVector> fieldValues;
0036 float Bx, By, Bz;
0037 int nLines = n1 * n2 * n3;
0038 fieldValues.reserve(nLines);
0039 for (int iLine = 0; iLine < nLines; ++iLine) {
0040 inFile >> Bx >> By >> Bz;
0041
0042
0043
0044
0045
0046
0047 Vector3DBase<double, LocalTag> lB = frame().toLocal(Vector3DBase<double, GlobalTag>(Bx, By, Bz));
0048 fieldValues.push_back(BVector(lB.x(), lB.y(), lB.z()));
0049 }
0050
0051 string lastEntry;
0052 inFile >> lastEntry;
0053 if (lastEntry != "complete") {
0054 cout << "ERROR during file reading: file is not complete" << endl;
0055 }
0056
0057 GlobalPoint grefp(GlobalPoint::Cylindrical(xref, yref, zref));
0058
0059 #ifdef DEBUG_GRID
0060 LocalPoint lrefp = frame().toLocal(grefp);
0061 cout << "Grid reference point in grid system: " << xref << "," << yref << "," << zref << endl;
0062 cout << "Grid reference point in global x,y,z: " << grefp << endl;
0063 cout << "Grid reference point in local x,y,z: " << lrefp << endl;
0064 cout << "steps " << stepx << "," << stepy << "," << stepz << endl;
0065 cout << "RParAsFunOfPhi[0...4] = ";
0066 for (int i = 0; i < 4; ++i)
0067 cout << RParAsFunOfPhi[i] << " ";
0068 cout << endl;
0069 #endif
0070
0071 Grid1D gridX(0, n1 - 1, n1);
0072 Grid1D gridY(yref, yref + stepy * (n2 - 1), n2);
0073 Grid1D gridZ(grefp.z(), grefp.z() + stepz * (n3 - 1), n3);
0074
0075 grid_ = GridType(gridX, gridY, gridZ, fieldValues);
0076
0077 stepConstTerm_ = (RParAsFunOfPhi[0] - RParAsFunOfPhi[2]) / (n1 - 1);
0078 stepPhiTerm_ = (RParAsFunOfPhi[1] - RParAsFunOfPhi[3]) / (n1 - 1);
0079 startConstTerm_ = RParAsFunOfPhi[2];
0080 startPhiTerm_ = RParAsFunOfPhi[3];
0081
0082
0083
0084
0085 }
0086
0087 MFGrid::LocalVector SpecialCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0088
0089
0090
0091 LinearGridInterpolator3D interpol(grid_);
0092 double a, b, c;
0093 toGridFrame(p, a, b, c);
0094
0095
0096
0097 return LocalVector(interpol.interpolate(a, b, c));
0098 }
0099
0100 void SpecialCylindricalMFGrid::dump() const {}
0101
0102 MFGrid::LocalPoint SpecialCylindricalMFGrid::fromGridFrame(double a, double b, double c) const {
0103 double sinPhi;
0104 if (sector1) {
0105 sinPhi = cos(b);
0106 } else {
0107 sinPhi = sin(b);
0108 }
0109
0110 double R = a * stepSize(sinPhi) + startingPoint(sinPhi);
0111
0112
0113 GlobalPoint gp(GlobalPoint::Cylindrical(R, b, c));
0114 return frame().toLocal(gp);
0115 }
0116
0117 void SpecialCylindricalMFGrid::toGridFrame(const LocalPoint& p, double& a, double& b, double& c) const {
0118 GlobalPoint gp = frame().toGlobal(p);
0119 double sinPhi;
0120 if (sector1) {
0121 sinPhi = cos(gp.phi());
0122 } else {
0123 sinPhi = sin(gp.phi());
0124 }
0125 a = (gp.perp() - startingPoint(sinPhi)) / stepSize(sinPhi);
0126
0127
0128 b = gp.phi();
0129 c = gp.z();
0130
0131 #ifdef DEBUG_GRID
0132 if (sector1) {
0133 cout << "toGridFrame: sinPhi ";
0134 } else {
0135 cout << "toGridFrame: cosPhi ";
0136 }
0137 cout << sinPhi << " LocalPoint " << p << " GlobalPoint " << gp << endl
0138 << " a " << a << " b " << b << " c " << c << endl;
0139
0140 #endif
0141 }