File indexing completed on 2021-02-14 13:31:24
0001 #include "prepareMagneticFieldGrid.h"
0002
0003 #include "MagneticField/Interpolation/src/VectorFieldInterpolation.h"
0004 #include "MagneticField/Interpolation/src/binary_ofstream.h"
0005 #include "DataFormats/GeometryVector/interface/Pi.h"
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <sstream>
0009 #include <iostream>
0010
0011 using namespace std;
0012
0013 void prepareMagneticFieldGrid::countTrueNumberOfPoints(const std::string &name) const {
0014 int nLines = 0;
0015 int nPoints = 0;
0016
0017
0018 IndexedDoubleVector XBVector;
0019 std::vector<IndexedDoubleVector> XBValues;
0020
0021
0022 double epsilonRadius = 1e-8;
0023 std::ifstream file(name.c_str());
0024 string line;
0025 if (file.good()) {
0026 while (getline(file, line)) {
0027 double x1, x2, x3, Bx, By, Bz, perm;
0028 stringstream linestr;
0029 linestr << line;
0030 linestr >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm;
0031 ;
0032 if (file) {
0033 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0034
0035
0036
0037
0038
0039 XBValues.push_back(XBVector);
0040 ++nLines;
0041 }
0042 }
0043 }
0044
0045
0046 for (int iLine = 0; iLine < nLines; ++iLine) {
0047 double x1, x2, x3, Bx, By, Bz;
0048 XBVector = XBValues.operator[](iLine);
0049 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0050 double pnt[3] = {x1, x2, x3};
0051 double tmpBz = Bz;
0052 bool isSinglePoint = true;
0053 for (int i = 0; i < nLines; ++i) {
0054 if (i < iLine) {
0055 XBVector = XBValues.operator[](i);
0056 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0057 double distPP =
0058 sqrt((pnt[0] - x1) * (pnt[0] - x1) + (pnt[1] - x2) * (pnt[1] - x2) + (pnt[2] - x3) * (pnt[2] - x3));
0059 if (distPP < epsilonRadius) {
0060 isSinglePoint = false;
0061 cout << "same points(" << iLine << ") with dR = " << distPP << endl;
0062 cout << " point 1: " << iLine << " " << pnt[0] << " " << pnt[1] << " " << pnt[2]
0063 << " Bz: " << tmpBz << endl;
0064 cout << " point 2: " << i << " " << x1 << " " << x2 << " " << x3 << " Bz: " << Bz << endl;
0065 }
0066 }
0067 }
0068 if (isSinglePoint)
0069 ++nPoints;
0070 }
0071
0072 if (PRINT)
0073 cout << " file " << name << endl;
0074 if (PRINT)
0075 cout << " # lines = " << nLines << " # points = " << nPoints;
0076 if (nLines == nPoints) {
0077 if (PRINT)
0078 cout << " --> PASSED" << endl;
0079 } else {
0080 if (PRINT)
0081 cout << " --> FAILED" << endl;
0082 }
0083
0084 return;
0085 }
0086
0087 void prepareMagneticFieldGrid::fillFromFile(const std::string &name) {
0088 double phi = (rotateSector - 1) * Geom::pi() / 6.;
0089 double sphi = sin(-phi);
0090 double cphi = cos(-phi);
0091
0092
0093 std::string::size_type ibeg, iend;
0094 ibeg = name.find('-');
0095 iend = name.rfind('-');
0096 if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "xyz") {
0097 XyzCoordinates = true;
0098 } else if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "rpz") {
0099 RpzCoordinates = true;
0100 } else {
0101 cout << " Unrecognized input file name: valid formats are *-xyz-* or *-rpz-*" << endl;
0102 abort();
0103 }
0104
0105
0106 IndexedDoubleVector XBVector;
0107 std::vector<IndexedDoubleVector> XBValues;
0108 std::vector<IndexedDoubleVector> IXBValues;
0109 std::vector<IndexedDoubleVector> XBArray;
0110
0111
0112 std::vector<double> valX[3];
0113 std::vector<int> numX[3];
0114
0115
0116 int nLines = 0;
0117 int index[3] = {0, 0, 0};
0118 int nSteps[3] = {0, 0, 0};
0119 double x1, x2, x3, Bx, By, Bz, perm, poten;
0120 std::ifstream file(name.c_str());
0121 int n_air = 0;
0122 int n_iron = 0;
0123
0124 if (file.good()) {
0125 while (file.good()) {
0126 file >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm >> poten;
0127 if (file) {
0128 if (rotateSector > 1) {
0129 if (RpzCoordinates) {
0130 x2 = x2 - phi;
0131 if (fabs(x2) > 0.78) {
0132 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << fabs(x2) << endl;
0133 abort();
0134 }
0135 } else {
0136 double x = x1;
0137 double y = x2;
0138 x1 = x * cphi - y * sphi;
0139 x2 = x * sphi + y * cphi;
0140 if (fabs(atan2(x2, x1)) > 0.78) {
0141 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << atan2(x2, x1)
0142 << endl;
0143 abort();
0144 }
0145 }
0146 double Bx0 = Bx;
0147 double By0 = By;
0148 Bx = Bx0 * cphi - By0 * sphi;
0149 By = Bx0 * sphi + By0 * cphi;
0150 }
0151
0152
0153
0154
0155 if (perm > 1.)
0156 ++n_iron;
0157 else
0158 ++n_air;
0159
0160 XBVector.putI3(index[0], index[1], index[2]);
0161 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0162 XBValues.push_back(XBVector);
0163 ++nLines;
0164
0165 double pnt[3] = {x1, x2, x3};
0166 if (nLines == 1) {
0167 for (int i = 0; i < 3; ++i) {
0168 valX[i].push_back(pnt[i]);
0169 numX[i].push_back(1);
0170 }
0171 } else {
0172 for (int i = 0; i < 3; ++i) {
0173 bool knownValue = false;
0174 for (int j = 0; j < (nSteps[i] + 1); ++j) {
0175 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0176 knownValue = true;
0177 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0178 numX[i].operator[](j) += 1;
0179 }
0180 if (!knownValue) {
0181 valX[i].push_back(pnt[i]);
0182 nSteps[i]++;
0183 numX[i].push_back(1);
0184 }
0185 }
0186 }
0187 }
0188 }
0189 if (n_air && n_iron)
0190 cout << "WARNING: table " << name << " contains " << n_air << " points with perm=1 and " << n_iron
0191 << " with perm>1" << endl;
0192 } else
0193 return;
0194
0195
0196
0197
0198 for (int i = 0; i < 3; ++i) {
0199
0200 sort(valX[i].begin(), valX[i].end());
0201
0202 EasyCoordinate[i] = true;
0203 NumberOfPoints[i] = nSteps[i] + 1;
0204 BasicDistance0[i] = 9e9;
0205 double newVal = 0.0;
0206 double oldVal = 0.0;
0207 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0208 newVal = valX[i].operator[](j);
0209 if (j == 1)
0210 BasicDistance0[i] = newVal - oldVal;
0211 else if (j > 1) {
0212 if (std::abs(newVal - oldVal - BasicDistance0[i]) > EPSILON)
0213 EasyCoordinate[i] = false;
0214 }
0215 oldVal = newVal;
0216 }
0217 }
0218
0219
0220 for (int iLine = 0; iLine < nLines; ++iLine) {
0221 XBVector = XBValues.operator[](iLine);
0222 XBVector.getI3(index[0], index[1], index[2]);
0223 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0224 double pnt[3] = {x1, x2, x3};
0225 for (int i = 0; i < 3; ++i) {
0226 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0227 if (EasyCoordinate[i]) {
0228 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0229 index[i] = j;
0230 }
0231 }
0232 }
0233 XBVector.putI3(index[0], index[1], index[2]);
0234 IXBValues.push_back(XBVector);
0235 }
0236 if (!XBValues.empty())
0237 XBValues.clear();
0238
0239
0240 sort(IXBValues.begin(), IXBValues.end());
0241
0242
0243 bool systematicGrid[3] = {false, false, false};
0244 KnownStructure = true;
0245 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines)
0246 KnownStructure = false;
0247 for (int i = 0; i < 3; ++i) {
0248 systematicGrid[i] = EasyCoordinate[i];
0249 if (!systematicGrid[i])
0250 KnownStructure = false;
0251 }
0252
0253 if (KnownStructure == true) {
0254 for (int iLine = 0; iLine < nLines; ++iLine) {
0255 XBVector = IXBValues.operator[](iLine);
0256 XBArray.push_back(XBVector);
0257 }
0258 }
0259
0260
0261 if (!XBArray.empty())
0262 XBVector = XBArray.front();
0263 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0264 double firstListPoint[3] = {x1, x2, x3};
0265 if (!XBArray.empty())
0266 XBVector = XBArray.back();
0267 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0268 double secondRefPoint[3] = {x1, x2, x3};
0269
0270
0271 KnownStructure = true;
0272 for (int i = 0; i < 3; ++i) {
0273 ReferencePoint[i] = firstListPoint[i];
0274
0275 double stepSize = BasicDistance0[i];
0276 for (int j = 0; j < 3; ++j) {
0277 stepSize += BasicDistance1[i][j] * nSteps[j];
0278 }
0279 double offset = 0.0;
0280 for (int j = 0; j < 3; ++j) {
0281 offset += BasicDistance2[i][j] * nSteps[j];
0282 }
0283
0284 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + offset + stepSize * nSteps[i]);
0285 if (std::abs(totDiff) < EPSILON * 10.)
0286 systematicGrid[i] = true;
0287 else
0288 systematicGrid[i] = false;
0289 if (!systematicGrid[i])
0290 KnownStructure = false;
0291 }
0292
0293 if (KnownStructure) {
0294 if (XyzCoordinates)
0295 GridType = 1;
0296 if (RpzCoordinates)
0297 GridType = 3;
0298 if (PRINT) {
0299
0300 cout << " read " << nLines << " lines -> grid structure:" << endl;
0301 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
0302 << " N3 = " << NumberOfPoints[2] << endl;
0303 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
0304 << " X3 = " << ReferencePoint[2] << endl;
0305 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
0306 << " A3 = " << BasicDistance0[2] << endl;
0307 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
0308 << " E3 = " << EasyCoordinate[2] << endl;
0309 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
0310 << " S3 = " << systematicGrid[2];
0311 if (KnownStructure)
0312 cout << " --> VALID (step I)" << endl;
0313 else
0314 cout << " --> INVALID (step I)" << endl;
0315 }
0316 }
0317
0318
0319
0320
0321 bool goodIndex[3] = {true, true, true};
0322
0323 if (!KnownStructure) {
0324
0325 int nEasy = 0;
0326 for (int i = 0; i < 3; ++i) {
0327 if (EasyCoordinate[i])
0328 ++nEasy;
0329 }
0330 std::vector<double> misValX[3];
0331
0332
0333 if (nEasy > 0) {
0334 bool lastLine = false;
0335 int lowerLineLimit = 0;
0336 int upperLineLimit = 0;
0337 int oldIndex[3] = {0, 0, 0};
0338 for (int iLine = 0; iLine < nLines; ++iLine) {
0339 if (iLine == (nLines - 1))
0340 lastLine = true;
0341 XBVector = IXBValues.operator[](iLine);
0342 XBVector.getI3(index[0], index[1], index[2]);
0343 bool newIndices = false;
0344 for (int i = 0; i < 3; ++i) {
0345 if (oldIndex[i] != index[i]) {
0346 oldIndex[i] = index[i];
0347 newIndices = true;
0348 }
0349 }
0350 if (newIndices) {
0351 for (int i = 0; i < 3; ++i) {
0352 if (!EasyCoordinate[i])
0353 sort(misValX[i].begin(), misValX[i].end());
0354 }
0355 lowerLineLimit = upperLineLimit;
0356 upperLineLimit = iLine;
0357 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0358 XBVector = IXBValues.operator[](jLine);
0359 XBVector.getI3(index[0], index[1], index[2]);
0360 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0361 double pnt[3] = {x1, x2, x3};
0362 for (int i = 0; i < 3; ++i) {
0363 if (!EasyCoordinate[i]) {
0364 double thisVal = 0.;
0365 double lastVal = 0.;
0366 double deltVal = 0.;
0367 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0368 thisVal = misValX[i].operator[](j);
0369 if (std::abs(thisVal - pnt[i]) < EPSILON)
0370 index[i] = j;
0371 if (j == 1)
0372 deltVal = thisVal - lastVal;
0373 else if (j > 1) {
0374 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0375 goodIndex[i] = false;
0376 }
0377 lastVal = thisVal;
0378 }
0379 }
0380 }
0381 XBVector.putI3(index[0], index[1], index[2]);
0382 XBArray.push_back(XBVector);
0383 for (int i = 0; i < 3; ++i) {
0384 if (!misValX[i].empty())
0385 misValX[i].clear();
0386 }
0387 }
0388 }
0389 XBVector = IXBValues.operator[](iLine);
0390 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0391 double pnt[3] = {x1, x2, x3};
0392 for (int i = 0; i < 3; ++i) {
0393 if (!EasyCoordinate[i]) {
0394 if (iLine == 0) {
0395 NumberOfPoints[i] = 1;
0396 misValX[i].push_back(pnt[i]);
0397 } else if (newIndices) {
0398 NumberOfPoints[i] = 1;
0399 misValX[i].push_back(pnt[i]);
0400 } else {
0401 bool knownValue = false;
0402 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0403 if (std::abs(misValX[i].operator[](j) - pnt[i]) < EPSILON)
0404 knownValue = true;
0405 }
0406 if (!knownValue) {
0407 ++NumberOfPoints[i];
0408 misValX[i].push_back(pnt[i]);
0409 }
0410 }
0411 }
0412 }
0413 if (lastLine) {
0414 for (int i = 0; i < 3; ++i) {
0415 if (!EasyCoordinate[i])
0416 sort(misValX[i].begin(), misValX[i].end());
0417 }
0418 lowerLineLimit = upperLineLimit;
0419 upperLineLimit = nLines;
0420 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0421 XBVector = IXBValues.operator[](jLine);
0422 XBVector.getI3(index[0], index[1], index[2]);
0423 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0424 double pnt[3] = {x1, x2, x3};
0425 for (int i = 0; i < 3; ++i) {
0426 if (!EasyCoordinate[i]) {
0427 double thisVal = 0.;
0428 double lastVal = 0.;
0429 double deltVal = 0.;
0430 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0431 thisVal = misValX[i].operator[](j);
0432 if (std::abs(thisVal - pnt[i]) < EPSILON)
0433 index[i] = j;
0434 if (j == 1)
0435 deltVal = thisVal - lastVal;
0436 else if (j > 1) {
0437 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0438 goodIndex[i] = false;
0439 }
0440 lastVal = thisVal;
0441 }
0442 }
0443 }
0444 XBVector.putI3(index[0], index[1], index[2]);
0445 XBArray.push_back(XBVector);
0446 for (int i = 0; i < 3; ++i) {
0447 if (!misValX[i].empty())
0448 misValX[i].clear();
0449 }
0450 }
0451 }
0452 }
0453 }
0454 int nGoodInd = 0;
0455 for (int i = 0; i < 3; ++i) {
0456 if (goodIndex[i])
0457 ++nGoodInd;
0458 }
0459
0460 if (nGoodInd < 3) {
0461 cout << endl;
0462 cout << "Neasy = 1: beginning of T E S T area!" << endl;
0463 cout << "# good indices = " << nGoodInd << " (" << goodIndex[0] << goodIndex[1] << goodIndex[2] << ") " << endl;
0464
0465 for (int iLine = 0; iLine < nLines; ++iLine) {
0466 XBVector = XBArray.operator[](iLine);
0467 XBVector.getI3(index[0], index[1], index[2]);
0468 for (int i = 0; i < 3; ++i) {
0469 if (!goodIndex[i])
0470 index[i] = 0;
0471 }
0472 XBVector.putI3(index[0], index[1], index[2]);
0473 }
0474 sort(XBArray.begin(), XBArray.end());
0475
0476 cout << "Neasy = 1: end of T E S T area!" << endl;
0477 cout << endl;
0478 }
0479
0480 if (nEasy > 0) {
0481
0482 sort(XBArray.begin(), XBArray.end());
0483
0484 double miniGrid[3][2][2][2] = {{{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}},
0485 {{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}},
0486 {{{0., 0.}, {0., 0.}}, {{0., 0.}, {0., 0.}}}};
0487 for (int iLine = 0; iLine < nLines; ++iLine) {
0488 XBVector = XBArray.operator[](iLine);
0489 XBVector.getI3(index[0], index[1], index[2]);
0490 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0491 double pnt[3] = {x1, x2, x3};
0492 if (index[0] < 2 && index[1] < 2 && index[2] < 2) {
0493 for (int i = 0; i < 3; ++i) {
0494 miniGrid[i][index[0]][index[1]][index[2]] = pnt[i];
0495 }
0496 }
0497 }
0498
0499
0500 BasicDistance0[0] = miniGrid[0][1][0][0] - miniGrid[0][0][0][0];
0501 BasicDistance0[1] = miniGrid[1][0][1][0] - miniGrid[1][0][0][0];
0502 BasicDistance0[2] = miniGrid[2][0][0][1] - miniGrid[2][0][0][0];
0503
0504 double disd10 = miniGrid[0][1][1][0] - miniGrid[0][0][1][0];
0505 double disd00 = BasicDistance0[0];
0506 double disd01 = miniGrid[0][1][0][1] - miniGrid[0][0][0][1];
0507 BasicDistance1[0][0] = 0.0;
0508 BasicDistance1[0][1] = disd10 - disd00;
0509 BasicDistance1[0][2] = disd01 - disd00;
0510 double dis1d0 = miniGrid[1][1][1][0] - miniGrid[1][1][0][0];
0511 double dis0d0 = BasicDistance0[1];
0512 double dis0d1 = miniGrid[1][0][1][1] - miniGrid[1][0][0][1];
0513 BasicDistance1[1][0] = dis1d0 - dis0d0;
0514 BasicDistance1[1][1] = 0.0;
0515 BasicDistance1[1][2] = dis0d1 - dis0d0;
0516 double dis10d = miniGrid[2][1][0][1] - miniGrid[2][1][0][0];
0517 double dis00d = BasicDistance0[2];
0518 double dis01d = miniGrid[2][0][1][1] - miniGrid[2][0][1][0];
0519 BasicDistance1[2][0] = dis10d - dis00d;
0520 BasicDistance1[2][1] = dis01d - dis00d;
0521 BasicDistance1[2][2] = 0.0;
0522
0523 BasicDistance2[0][0] = 0.0;
0524 BasicDistance2[0][1] = miniGrid[0][0][1][0] - miniGrid[0][0][0][0];
0525 BasicDistance2[0][2] = miniGrid[0][0][0][1] - miniGrid[0][0][0][0];
0526 BasicDistance2[1][0] = miniGrid[1][1][0][0] - miniGrid[1][0][0][0];
0527 BasicDistance2[1][1] = 0.0;
0528 BasicDistance2[1][2] = miniGrid[1][0][0][1] - miniGrid[1][0][0][0];
0529 BasicDistance2[2][0] = miniGrid[2][1][0][0] - miniGrid[2][0][0][0];
0530 BasicDistance2[2][1] = miniGrid[2][0][1][0] - miniGrid[2][0][0][0];
0531 BasicDistance2[2][2] = 0.0;
0532
0533 for (int i = 0; i < 3; ++i) {
0534 if (NumberOfPoints[i] < 2)
0535 BasicDistance0[i] = 9e9;
0536 }
0537 }
0538
0539
0540 if (!XBArray.empty())
0541 XBVector = XBArray.front();
0542 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0543 firstListPoint[0] = x1;
0544 firstListPoint[1] = x2;
0545 firstListPoint[2] = x3;
0546 if (!XBArray.empty())
0547 XBVector = XBArray.back();
0548 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0549 secondRefPoint[0] = x1;
0550 secondRefPoint[1] = x2;
0551 secondRefPoint[2] = x3;
0552
0553
0554 KnownStructure = true;
0555 for (int i = 0; i < 3; ++i) {
0556 ReferencePoint[i] = firstListPoint[i];
0557 nSteps[i] = NumberOfPoints[i] - 1;
0558
0559 double stepSize = BasicDistance0[i];
0560 for (int j = 0; j < 3; ++j) {
0561 stepSize += BasicDistance1[i][j] * nSteps[j];
0562 }
0563 double offset = 0.0;
0564 for (int j = 0; j < 3; ++j) {
0565 offset += BasicDistance2[i][j] * nSteps[j];
0566 }
0567
0568 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + offset + stepSize * nSteps[i]);
0569 if (std::abs(totDiff) < EPSILON * 10.)
0570 systematicGrid[i] = true;
0571 else
0572 systematicGrid[i] = false;
0573 if (!systematicGrid[i])
0574 KnownStructure = false;
0575 }
0576
0577 if (KnownStructure) {
0578 if (XyzCoordinates)
0579 GridType = 2;
0580 if (RpzCoordinates)
0581 GridType = 4;
0582 if (PRINT) {
0583
0584 cout << " read " << nLines << " lines -> grid structure:" << endl;
0585 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
0586 << " N3 = " << NumberOfPoints[2] << endl;
0587 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
0588 << " X3 = " << ReferencePoint[2] << endl;
0589 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
0590 << " A3 = " << BasicDistance0[2] << endl;
0591 cout << " step parm1 : B11 = " << BasicDistance1[0][0] << " B21 = " << BasicDistance1[1][0]
0592 << " B31 = " << BasicDistance1[2][0] << endl;
0593 cout << " : B12 = " << BasicDistance1[0][1] << " B22 = " << BasicDistance1[1][1]
0594 << " B32 = " << BasicDistance1[2][1] << endl;
0595 cout << " : B13 = " << BasicDistance1[0][2] << " B23 = " << BasicDistance1[1][2]
0596 << " B33 = " << BasicDistance1[2][2] << endl;
0597 cout << " offset parm: O11 = " << BasicDistance2[0][0] << " O21 = " << BasicDistance2[1][0]
0598 << " O31 = " << BasicDistance2[2][0] << endl;
0599 cout << " : O12 = " << BasicDistance2[0][1] << " O22 = " << BasicDistance2[1][1]
0600 << " O32 = " << BasicDistance2[2][1] << endl;
0601 cout << " : O13 = " << BasicDistance2[0][2] << " O23 = " << BasicDistance2[1][2]
0602 << " O33 = " << BasicDistance2[2][2] << endl;
0603 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
0604 << " E3 = " << EasyCoordinate[2] << endl;
0605 cout << " : I1 = " << goodIndex[0] << " I2 = " << goodIndex[1] << " I3 = " << goodIndex[2]
0606 << endl;
0607 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
0608 << " S3 = " << systematicGrid[2];
0609 if (KnownStructure)
0610 cout << " --> VALID (step II)" << endl;
0611 else {
0612 cout << " --> INVALID (step II)" << endl;
0613 cout << " reason for error: ";
0614 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines) {
0615 cout << endl;
0616 cout << " N1*N2*N3 =/= N.lines --> exiting now ..." << endl;
0617 return;
0618 } else {
0619 cout << " no idea so far --> exiting now ..." << endl;
0620 return;
0621 }
0622 }
0623 }
0624 }
0625 }
0626
0627
0628
0629 SixDPoint GridEntry;
0630 for (unsigned int iLine = 0; iLine < XBArray.size(); ++iLine) {
0631 XBVector = XBArray.operator[](iLine);
0632 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0633 GridEntry.putP6(x1, x2, x3, Bx, By, Bz);
0634 GridData.push_back(GridEntry);
0635 }
0636 if (!XBArray.empty())
0637 XBArray.clear();
0638
0639 return;
0640 }
0641
0642 void prepareMagneticFieldGrid::fillFromFileSpecial(const std::string &name) {
0643
0644 GridType = 6;
0645
0646 double phi = (rotateSector - 1) * Geom::pi() / 6.;
0647 double sphi = sin(-phi);
0648 double cphi = cos(-phi);
0649
0650
0651 std::string::size_type ibeg, iend;
0652 ibeg = name.find('-');
0653 iend = name.rfind('-');
0654 if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "xyz") {
0655 XyzCoordinates = true;
0656 } else if ((name.substr(ibeg + 1, iend - ibeg - 1)) == "rpz") {
0657 RpzCoordinates = true;
0658 } else {
0659 cout << " Unrecognized input file name: valid formats are *-xyz-* or *-rpz-*" << endl;
0660 abort();
0661 }
0662
0663
0664 IndexedDoubleVector XBVector;
0665 std::vector<IndexedDoubleVector> XBValues;
0666 std::vector<IndexedDoubleVector> IXBValues;
0667 std::vector<IndexedDoubleVector> XBArray;
0668
0669
0670 std::vector<double> valX[3];
0671 std::vector<int> numX[3];
0672
0673
0674 int nLines = 0;
0675 int index[3] = {0, 0, 0};
0676 int nSteps[3] = {0, 0, 0};
0677 double x1, x2, x3, Bx, By, Bz, perm, poten;
0678 std::ifstream file(name.c_str());
0679 if (file.good()) {
0680 while (file.good()) {
0681 file >> x1 >> x2 >> x3 >> Bx >> By >> Bz >> perm >> poten;
0682 if (file) {
0683 if (rotateSector > 1) {
0684 if (RpzCoordinates) {
0685 x2 = x2 - phi;
0686 if (fabs(x2) > 0.78) {
0687 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << x2 << endl;
0688 abort();
0689 }
0690 } else {
0691 double x = x1;
0692 double y = x2;
0693 x1 = x * cphi - y * sphi;
0694 x2 = x * sphi + y * cphi;
0695 if (fabs(atan2(x2, x1)) > 0.78) {
0696 cout << "ERROR: Input coordinates do not belong to sector " << rotateSector << " : " << atan2(x2, x1)
0697 << endl;
0698 abort();
0699 }
0700 }
0701 double Bx0 = Bx;
0702 double By0 = By;
0703 Bx = Bx0 * cphi - By0 * sphi;
0704 By = Bx0 * sphi + By0 * cphi;
0705 }
0706
0707 XBVector.putI3(index[0], index[1], index[2]);
0708 XBVector.putV6(x1, x2, x3, Bx, By, Bz);
0709 XBValues.push_back(XBVector);
0710 ++nLines;
0711
0712
0713 double pnt[3] = {x1, x2, x3};
0714 if (nLines == 1) {
0715 for (int i = 0; i < 3; ++i) {
0716 valX[i].push_back(pnt[i]);
0717 numX[i].push_back(1);
0718 }
0719 } else {
0720 for (int i = 0; i < 3; ++i) {
0721 bool knownValue = false;
0722 for (int j = 0; j < (nSteps[i] + 1); ++j) {
0723 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0724 knownValue = true;
0725 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0726 numX[i].operator[](j) += 1;
0727 }
0728 if (!knownValue) {
0729 valX[i].push_back(pnt[i]);
0730 nSteps[i]++;
0731 numX[i].push_back(1);
0732 }
0733 }
0734 }
0735 }
0736 }
0737 } else
0738 return;
0739
0740
0741 for (int i = 0; i < 3; ++i) {
0742
0743 sort(valX[i].begin(), valX[i].end());
0744
0745 EasyCoordinate[i] = true;
0746 NumberOfPoints[i] = nSteps[i] + 1;
0747 BasicDistance0[i] = 9e9;
0748 double newVal = 0.0;
0749 double oldVal = 0.0;
0750 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0751 newVal = valX[i].operator[](j);
0752 if (j == 1)
0753 BasicDistance0[i] = newVal - oldVal;
0754 else if (j > 1) {
0755 if (std::abs(newVal - oldVal - BasicDistance0[i]) > EPSILON)
0756 EasyCoordinate[i] = false;
0757 }
0758 oldVal = newVal;
0759 }
0760 }
0761
0762 std::cout << std::endl;
0763
0764
0765 for (int iLine = 0; iLine < nLines; ++iLine) {
0766 XBVector = XBValues.operator[](iLine);
0767 XBVector.getI3(index[0], index[1], index[2]);
0768 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0769 double pnt[3] = {x1, x2, x3};
0770 for (int i = 0; i < 3; ++i) {
0771 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0772 if (EasyCoordinate[i]) {
0773 if (std::abs(valX[i].operator[](j) - pnt[i]) < EPSILON)
0774 index[i] = j;
0775 }
0776 }
0777 }
0778 XBVector.putI3(index[0], index[1], index[2]);
0779 IXBValues.push_back(XBVector);
0780 }
0781 if (!XBValues.empty())
0782 XBValues.clear();
0783
0784
0785 sort(IXBValues.begin(), IXBValues.end());
0786
0787
0788 bool systematicGrid[3] = {false, false, false};
0789 KnownStructure = true;
0790 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines)
0791 KnownStructure = false;
0792 for (int i = 0; i < 3; ++i) {
0793 systematicGrid[i] = EasyCoordinate[i];
0794 if (!systematicGrid[i])
0795 KnownStructure = false;
0796 }
0797
0798 if (KnownStructure == true) {
0799 for (int iLine = 0; iLine < nLines; ++iLine) {
0800 XBVector = IXBValues.operator[](iLine);
0801 XBArray.push_back(XBVector);
0802 }
0803 }
0804
0805
0806 if (!XBArray.empty())
0807 XBVector = XBArray.front();
0808 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0809 double firstListPoint[3] = {x1, x2, x3};
0810 if (!XBArray.empty())
0811 XBVector = XBArray.back();
0812 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0813 double secondRefPoint[3] = {x1, x2, x3};
0814
0815
0816 KnownStructure = true;
0817 for (int i = 0; i < 3; ++i) {
0818 ReferencePoint[i] = firstListPoint[i];
0819
0820 double stepSize = BasicDistance0[i];
0821
0822 double totDiff = secondRefPoint[i] - (ReferencePoint[i] + stepSize * nSteps[i]);
0823 if (std::abs(totDiff) < EPSILON * 10.)
0824 systematicGrid[i] = true;
0825 else
0826 systematicGrid[i] = false;
0827 if (!systematicGrid[i])
0828 KnownStructure = false;
0829 }
0830 bool goodIndex[3] = {true, true, true};
0831
0832 if (!KnownStructure) {
0833
0834 int nEasy = 0;
0835 for (int i = 0; i < 3; ++i) {
0836 if (EasyCoordinate[i])
0837 ++nEasy;
0838 }
0839 std::vector<double> misValX[3];
0840
0841
0842 if (nEasy > 0) {
0843 bool lastLine = false;
0844 int lowerLineLimit = 0;
0845 int upperLineLimit = 0;
0846 int oldIndex[3] = {0, 0, 0};
0847 for (int iLine = 0; iLine < nLines; ++iLine) {
0848 if (iLine == (nLines - 1))
0849 lastLine = true;
0850 XBVector = IXBValues.operator[](iLine);
0851 XBVector.getI3(index[0], index[1], index[2]);
0852 bool newIndices = false;
0853 for (int i = 0; i < 3; ++i) {
0854 if (oldIndex[i] != index[i]) {
0855 oldIndex[i] = index[i];
0856 newIndices = true;
0857 }
0858 }
0859 if (newIndices) {
0860 for (int i = 0; i < 3; ++i) {
0861 if (!EasyCoordinate[i])
0862 sort(misValX[i].begin(), misValX[i].end());
0863 }
0864 lowerLineLimit = upperLineLimit;
0865 upperLineLimit = iLine;
0866 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0867 XBVector = IXBValues.operator[](jLine);
0868 XBVector.getI3(index[0], index[1], index[2]);
0869 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0870 double pnt[3] = {x1, x2, x3};
0871 for (int i = 0; i < 3; ++i) {
0872 if (!EasyCoordinate[i]) {
0873 double thisVal = 0.;
0874 double lastVal = 0.;
0875 double deltVal = 0.;
0876 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0877 thisVal = misValX[i].operator[](j);
0878 if (std::abs(thisVal - pnt[i]) < EPSILON)
0879 index[i] = j;
0880 if (j == 1)
0881 deltVal = thisVal - lastVal;
0882 else if (j > 1) {
0883 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0884 goodIndex[i] = false;
0885 }
0886 lastVal = thisVal;
0887 }
0888 }
0889 }
0890 XBVector.putI3(index[0], index[1], index[2]);
0891 XBArray.push_back(XBVector);
0892 for (int i = 0; i < 3; ++i) {
0893 if (!misValX[i].empty())
0894 misValX[i].clear();
0895 }
0896 }
0897 }
0898 XBVector = IXBValues.operator[](iLine);
0899 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0900 double pnt[3] = {x1, x2, x3};
0901 for (int i = 0; i < 3; ++i) {
0902 if (!EasyCoordinate[i]) {
0903 if (iLine == 0) {
0904 NumberOfPoints[i] = 1;
0905 misValX[i].push_back(pnt[i]);
0906 } else if (newIndices) {
0907 NumberOfPoints[i] = 1;
0908 misValX[i].push_back(pnt[i]);
0909 } else {
0910 bool knownValue = false;
0911 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0912 if (std::abs(misValX[i].operator[](j) - pnt[i]) < EPSILON)
0913 knownValue = true;
0914 }
0915 if (!knownValue) {
0916 ++NumberOfPoints[i];
0917 misValX[i].push_back(pnt[i]);
0918 }
0919 }
0920 }
0921 }
0922 if (lastLine) {
0923 for (int i = 0; i < 3; ++i) {
0924 if (!EasyCoordinate[i])
0925 sort(misValX[i].begin(), misValX[i].end());
0926 }
0927 lowerLineLimit = upperLineLimit;
0928 upperLineLimit = nLines;
0929 for (int jLine = lowerLineLimit; jLine < upperLineLimit; ++jLine) {
0930 XBVector = IXBValues.operator[](jLine);
0931 XBVector.getI3(index[0], index[1], index[2]);
0932 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0933 double pnt[3] = {x1, x2, x3};
0934 for (int i = 0; i < 3; ++i) {
0935 if (!EasyCoordinate[i]) {
0936 double thisVal = 0.;
0937 double lastVal = 0.;
0938 double deltVal = 0.;
0939 for (int j = 0; j < NumberOfPoints[i]; ++j) {
0940 thisVal = misValX[i].operator[](j);
0941 if (std::abs(thisVal - pnt[i]) < EPSILON)
0942 index[i] = j;
0943 if (j == 1)
0944 deltVal = thisVal - lastVal;
0945 else if (j > 1) {
0946 if (std::abs(thisVal - lastVal - deltVal) > EPSILON)
0947 goodIndex[i] = false;
0948 }
0949 lastVal = thisVal;
0950 }
0951 }
0952 }
0953 XBVector.putI3(index[0], index[1], index[2]);
0954 XBArray.push_back(XBVector);
0955 for (int i = 0; i < 3; ++i) {
0956 if (!misValX[i].empty())
0957 misValX[i].clear();
0958 }
0959 }
0960 }
0961 }
0962 }
0963 int nGoodInd = 0;
0964 for (int i = 0; i < 3; ++i) {
0965 if (goodIndex[i])
0966 ++nGoodInd;
0967 }
0968
0969 if (nEasy > 0) {
0970
0971 sort(XBArray.begin(), XBArray.end());
0972
0973 for (int i = 0; i < 3; ++i) {
0974 nSteps[i] = NumberOfPoints[i] - 1;
0975 }
0976 std::vector<double> rMinVec;
0977 std::vector<double> rMaxVec;
0978 std::vector<double> phiVec;
0979 for (int iLine = 0; iLine < nLines; ++iLine) {
0980 XBVector = XBArray.operator[](iLine);
0981 XBVector.getI3(index[0], index[1], index[2]);
0982
0983 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
0984 if (index[2] == 0) {
0985 if (index[0] == 0)
0986 rMinVec.push_back(x1);
0987 if (index[0] == nSteps[0])
0988 rMaxVec.push_back(x1);
0989 if (index[0] == nSteps[0])
0990 phiVec.push_back(x2);
0991 }
0992 }
0993 for (int j = 0; j < NumberOfPoints[1]; ++j) {
0994 double phi = phiVec.operator[](j);
0995 double rMin = rMinVec.operator[](j);
0996 double rMax = rMaxVec.operator[](j);
0997 double rhoMin, rhoMax;
0998 if (GridType == 6) {
0999 rhoMin = rMin * cos(phi);
1000 rhoMax = rMax * cos(phi);
1001 } else if (GridType == 5) {
1002 rhoMin = rMin * sin(phi);
1003 rhoMax = rMax * sin(phi);
1004 } else {
1005 cout << "unknown grid type!" << endl;
1006 abort();
1007 }
1008
1009 if (j == 0) {
1010 RParAsFunOfPhi[0] = rMax;
1011 RParAsFunOfPhi[1] = rhoMax;
1012 RParAsFunOfPhi[2] = rMin;
1013 RParAsFunOfPhi[3] = rhoMin;
1014 } else {
1015 if (std::abs(rMax - RParAsFunOfPhi[0]) > EPSILON)
1016 RParAsFunOfPhi[0] = 0.;
1017 if (std::abs(rhoMax - RParAsFunOfPhi[1]) > EPSILON)
1018 RParAsFunOfPhi[1] = 0.;
1019 if (std::abs(rMin - RParAsFunOfPhi[2]) > EPSILON)
1020 RParAsFunOfPhi[2] = 0.;
1021 if (std::abs(rhoMin - RParAsFunOfPhi[3]) > EPSILON)
1022 RParAsFunOfPhi[3] = 0.;
1023 }
1024 }
1025
1026
1027 for (int i = 0; i < 3; ++i) {
1028 if (NumberOfPoints[i] < 2)
1029 BasicDistance0[i] = 9e9;
1030 }
1031 }
1032
1033
1034 if (!XBArray.empty())
1035 XBVector = XBArray.front();
1036 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1037 firstListPoint[0] = x1;
1038 firstListPoint[1] = x2;
1039 firstListPoint[2] = x3;
1040 if (!XBArray.empty())
1041 XBVector = XBArray.back();
1042 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1043 secondRefPoint[0] = x1;
1044 secondRefPoint[1] = x2;
1045 secondRefPoint[2] = x3;
1046
1047
1048 KnownStructure = true;
1049 for (int i = 0; i < 3; ++i) {
1050 ReferencePoint[i] = firstListPoint[i];
1051 nSteps[i] = NumberOfPoints[i] - 1;
1052 }
1053
1054 double sinPhi;
1055 if (GridType == 6) {
1056 sinPhi = cos(secondRefPoint[1]);
1057 } else if (GridType == 5) {
1058 sinPhi = sin(secondRefPoint[1]);
1059 } else {
1060 cout << "unknown GridType!" << endl;
1061 abort();
1062 }
1063
1064 if (std::abs(sinPhi) < EPSILON) {
1065 cout << " WARNING: unexpected sinPhi parameter = 0" << endl;
1066 sinPhi = EPSILON;
1067 }
1068
1069 double totStepSize =
1070 RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1071 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1072 double totDiff = secondRefPoint[0] - (startingPoint + totStepSize);
1073 if (std::abs(totDiff) < EPSILON * 10.)
1074 systematicGrid[0] = true;
1075 else
1076 systematicGrid[0] = false;
1077 if (!systematicGrid[0])
1078 KnownStructure = false;
1079 totDiff = secondRefPoint[1] - (ReferencePoint[1] + BasicDistance0[1] * nSteps[1]);
1080 if (std::abs(totDiff) < EPSILON * 10.)
1081 systematicGrid[1] = true;
1082 else
1083 systematicGrid[1] = false;
1084 if (!systematicGrid[1])
1085 KnownStructure = false;
1086 totDiff = secondRefPoint[2] - (ReferencePoint[2] + BasicDistance0[2] * nSteps[2]);
1087 if (std::abs(totDiff) < EPSILON * 10.)
1088 systematicGrid[2] = true;
1089 else
1090 systematicGrid[2] = false;
1091 if (!systematicGrid[2])
1092 KnownStructure = false;
1093
1094 if (PRINT) {
1095
1096 cout << " read " << nLines << " lines -> grid structure:" << endl;
1097 cout << " # of points: N1 = " << NumberOfPoints[0] << " N2 = " << NumberOfPoints[1]
1098 << " N3 = " << NumberOfPoints[2] << endl;
1099 cout << " ref. point : X1 = " << ReferencePoint[0] << " X2 = " << ReferencePoint[1]
1100 << " X3 = " << ReferencePoint[2] << endl;
1101 cout << " step parm0 : A1 = " << BasicDistance0[0] << " A2 = " << BasicDistance0[1]
1102 << " A3 = " << BasicDistance0[2] << endl;
1103 cout << " r param.s. : R1 = " << RParAsFunOfPhi[0] << " RHO1 = " << RParAsFunOfPhi[1]
1104 << " R2 = " << RParAsFunOfPhi[2] << " RHO2 = " << RParAsFunOfPhi[3] << endl;
1105 cout << " easy grid : E1 = " << EasyCoordinate[0] << " E2 = " << EasyCoordinate[1]
1106 << " E3 = " << EasyCoordinate[2] << endl;
1107 cout << " : I1 = " << goodIndex[0] << " I2 = " << goodIndex[1] << " I3 = " << goodIndex[2]
1108 << endl;
1109 cout << " structure : S1 = " << systematicGrid[0] << " S2 = " << systematicGrid[1]
1110 << " S3 = " << systematicGrid[2];
1111 if (KnownStructure)
1112 cout << " --> VALID (step II)" << endl;
1113 else {
1114 cout << " --> INVALID (step II)" << endl;
1115 cout << " reason for error: ";
1116 if (NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2] != nLines) {
1117 cout << endl;
1118 cout << " N1*N2*N3 =/= N.lines --> exiting now ..." << endl;
1119 return;
1120 } else {
1121 cout << " no idea so far --> exiting now ..." << endl;
1122 return;
1123 }
1124 }
1125 }
1126 }
1127
1128
1129 SixDPoint GridEntry;
1130 for (unsigned int iLine = 0; iLine < XBArray.size(); ++iLine) {
1131 XBVector = XBArray.operator[](iLine);
1132 XBVector.getV6(x1, x2, x3, Bx, By, Bz);
1133 GridEntry.putP6(x1, x2, x3, Bx, By, Bz);
1134 GridData.push_back(GridEntry);
1135 }
1136 if (!XBArray.empty())
1137 XBArray.clear();
1138
1139 return;
1140 }
1141
1142 int prepareMagneticFieldGrid::gridType() {
1143 int type = GridType;
1144 if (PRINT) {
1145 if (type == 0)
1146 cout << " grid type = " << type << " --> not determined" << endl;
1147 if (type == 1)
1148 cout << " grid type = " << type << " --> (x,y,z) cube" << endl;
1149 if (type == 2)
1150 cout << " grid type = " << type << " --> (x,y,z) trapezoid" << endl;
1151 if (type == 3)
1152 cout << " grid type = " << type << " --> (r,phi,z) cube" << endl;
1153 if (type == 4)
1154 cout << " grid type = " << type << " --> (r,phi,z) trapezoid" << endl;
1155 if (type == 5)
1156 cout << " grid type = " << type << " --> (r,phi,z) 1/sin(phi)" << endl;
1157 if (type == 6)
1158 cout << " grid type = " << type << " --> (r,phi,z) 1/cos(phi)" << endl;
1159 }
1160 return type;
1161 }
1162
1163 void prepareMagneticFieldGrid::validateAllPoints() {
1164 if (!KnownStructure)
1165 return;
1166
1167 double x1, x2, x3, Bx, By, Bz;
1168 int numberOfErrors = 0;
1169
1170
1171 int index[3];
1172
1173 if (GridType < 5) {
1174 for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1175 for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1176 for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1177
1178 putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1179
1180
1181 double pnt[3] = {x1, x2, x3};
1182 int tmpIdx[3] = {999, 999, 999};
1183 putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1184 for (int i = 0; i < 3; ++i) {
1185 if (tmpIdx[i] != index[i])
1186 ++numberOfErrors;
1187 }
1188
1189
1190 double tmpPnt[3] = {0.0, 0.0, 0.0};
1191 putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1192 for (int i = 0; i < 3; ++i) {
1193 if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1194 ++numberOfErrors;
1195 }
1196 }
1197 }
1198 }
1199 }
1200
1201 if (GridType == 5 || GridType == 6) {
1202 for (index[0] = 0; index[0] < NumberOfPoints[0]; ++index[0]) {
1203 for (index[1] = 0; index[1] < NumberOfPoints[1]; ++index[1]) {
1204 for (index[2] = 0; index[2] < NumberOfPoints[2]; ++index[2]) {
1205
1206 putIndicesGetXAndB(index[0], index[1], index[2], x1, x2, x3, Bx, By, Bz);
1207
1208
1209 double pnt[3] = {x1, x2, x3};
1210 int tmpIdx[3] = {999, 999, 999};
1211 putCoordGetIndices((x1 + EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1212 for (int i = 0; i < 3; ++i) {
1213 if (tmpIdx[i] != index[i]) {
1214 putCoordGetIndices((x1 - EPSILON), (x2 + EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1215 if (tmpIdx[i] != index[i]) {
1216 putCoordGetIndices((x1 + EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1217 if (tmpIdx[i] != index[i]) {
1218 putCoordGetIndices((x1 - EPSILON), (x2 - EPSILON), (x3 + EPSILON), tmpIdx[0], tmpIdx[1], tmpIdx[2]);
1219 if (tmpIdx[i] != index[i])
1220 ++numberOfErrors;
1221 }
1222 }
1223 }
1224 }
1225
1226
1227 double tmpPnt[3] = {0.0, 0.0, 0.0};
1228 putIndCalcXReturnB(index[0], index[1], index[2], tmpPnt[0], tmpPnt[1], tmpPnt[2], Bx, By, Bz);
1229 for (int i = 0; i < 3; ++i) {
1230 if (std::abs(tmpPnt[i] - pnt[i]) > EPSILON)
1231 ++numberOfErrors;
1232 }
1233 }
1234 }
1235 }
1236 }
1237
1238 if (numberOfErrors > 0)
1239 KnownStructure = false;
1240 if (PRINT)
1241 cout << " grid validation of all points done --> # errors = " << numberOfErrors << endl;
1242
1243 return;
1244 }
1245
1246 void prepareMagneticFieldGrid::saveGridToFile(const std::string &outName) {
1247
1248 binary_ofstream outFile(outName);
1249
1250 outFile << GridType;
1251
1252 convertUnits();
1253 if (GridType == 1) {
1254 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1255 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1256 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1257 }
1258 if (GridType == 2) {
1259 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1260 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1261 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1262 outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1263 outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1264 outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1265 outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1266 outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1267 outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1268 outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1269 }
1270 if (GridType == 3) {
1271 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1272 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1273 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1274 }
1275 if (GridType == 4) {
1276 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1277 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1278 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1279 outFile << BasicDistance1[0][0] << BasicDistance1[1][0] << BasicDistance1[2][0];
1280 outFile << BasicDistance1[0][1] << BasicDistance1[1][1] << BasicDistance1[2][1];
1281 outFile << BasicDistance1[0][2] << BasicDistance1[1][2] << BasicDistance1[2][2];
1282 outFile << BasicDistance2[0][0] << BasicDistance2[1][0] << BasicDistance2[2][0];
1283 outFile << BasicDistance2[0][1] << BasicDistance2[1][1] << BasicDistance2[2][1];
1284 outFile << BasicDistance2[0][2] << BasicDistance2[1][2] << BasicDistance2[2][2];
1285 outFile << EasyCoordinate[0] << EasyCoordinate[1] << EasyCoordinate[2];
1286 }
1287 if (GridType == 5 || GridType == 6) {
1288 outFile << NumberOfPoints[0] << NumberOfPoints[1] << NumberOfPoints[2];
1289 outFile << ReferencePoint[0] << ReferencePoint[1] << ReferencePoint[2];
1290 outFile << BasicDistance0[0] << BasicDistance0[1] << BasicDistance0[2];
1291 outFile << RParAsFunOfPhi[0] << RParAsFunOfPhi[1] << RParAsFunOfPhi[2] << RParAsFunOfPhi[3];
1292 }
1293 int nLines = NumberOfPoints[0] * NumberOfPoints[1] * NumberOfPoints[2];
1294 SixDPoint GridPoint;
1295
1296 for (int iLine = 0; iLine < nLines; ++iLine) {
1297 GridPoint = GridData.operator[](iLine);
1298 float Bx = float(GridPoint.bx());
1299 float By = float(GridPoint.by());
1300 float Bz = float(GridPoint.bz());
1301 outFile << Bx << By << Bz;
1302
1303 }
1304
1305 const std::string lastEntry = "complete";
1306 outFile << lastEntry;
1307 outFile.close();
1308 if (PRINT)
1309 cout << " output " << outName << endl;
1310 return;
1311 }
1312
1313 void prepareMagneticFieldGrid::convertUnits() {
1314 double cm = 100.;
1315 if (XyzCoordinates) {
1316 for (int i = 0; i < 3; ++i) {
1317 ReferencePoint[i] *= cm;
1318 };
1319 for (int i = 0; i < 3; ++i) {
1320 BasicDistance0[i] *= cm;
1321 };
1322 for (int i = 0; i < 3; ++i) {
1323 for (int j = 0; j < 3; ++j) {
1324 BasicDistance1[i][j] *= cm;
1325 };
1326 };
1327 for (int i = 0; i < 3; ++i) {
1328 for (int j = 0; j < 3; ++j) {
1329 BasicDistance2[i][j] *= cm;
1330 };
1331 };
1332 for (int i = 0; i < 4; ++i) {
1333 RParAsFunOfPhi[i] *= cm;
1334 };
1335 }
1336 double du[3] = {100., 1., 100.};
1337 if (RpzCoordinates) {
1338 for (int i = 0; i < 3; ++i) {
1339 ReferencePoint[i] *= du[i];
1340 };
1341 for (int i = 0; i < 3; ++i) {
1342 BasicDistance0[i] *= du[i];
1343 };
1344 for (int i = 0; i < 3; ++i) {
1345 for (int j = 0; j < 3; ++j) {
1346 BasicDistance1[i][j] *= du[i];
1347 };
1348 };
1349 for (int i = 0; i < 3; ++i) {
1350 for (int j = 0; j < 3; ++j) {
1351 BasicDistance2[i][j] *= du[i];
1352 };
1353 };
1354 for (int i = 0; i < 4; ++i) {
1355 RParAsFunOfPhi[i] *= cm;
1356 };
1357 }
1358 return;
1359 }
1360
1361 bool prepareMagneticFieldGrid::isReady() { return KnownStructure; }
1362
1363 void prepareMagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, double &Bx, double &By, double &Bz) {
1364 Bx = By = Bz = 0.0;
1365 if (KnownStructure) {
1366
1367 VectorFieldInterpolation MagInterpol;
1368
1369 int index[3];
1370 putCoordGetIndices(X1, X2, X3, index[0], index[1], index[2]);
1371 int index0[3] = {0, 0, 0};
1372 int index1[3] = {0, 0, 0};
1373 for (int i = 0; i < 3; ++i) {
1374 if (NumberOfPoints[i] > 1) {
1375 index0[i] = std::max(0, index[i]);
1376 if (index0[i] > NumberOfPoints[i] - 2)
1377 index0[i] = NumberOfPoints[i] - 2;
1378 ;
1379 index1[i] = std::max(1, index[i] + 1);
1380 ;
1381 if (index1[i] > NumberOfPoints[i] - 1)
1382 index1[i] = NumberOfPoints[i] - 1;
1383 }
1384 }
1385 double tmpX[3];
1386 double tmpB[3];
1387
1388 putIndicesGetXAndB(index0[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1389 MagInterpol.defineCellPoint000(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1390 putIndicesGetXAndB(index1[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1391 MagInterpol.defineCellPoint100(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1392 putIndicesGetXAndB(index0[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1393 MagInterpol.defineCellPoint010(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1394 putIndicesGetXAndB(index1[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1395 MagInterpol.defineCellPoint110(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1396 putIndicesGetXAndB(index0[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1397 MagInterpol.defineCellPoint001(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1398 putIndicesGetXAndB(index1[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1399 MagInterpol.defineCellPoint101(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1400 putIndicesGetXAndB(index0[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1401 MagInterpol.defineCellPoint011(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1402 putIndicesGetXAndB(index1[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1403 MagInterpol.defineCellPoint111(tmpX[0], tmpX[1], tmpX[2], tmpB[0], tmpB[1], tmpB[2]);
1404
1405 MagInterpol.putSCoordGetVField(X1, X2, X3, Bx, By, Bz);
1406 }
1407
1408 return;
1409 }
1410
1411 void prepareMagneticFieldGrid::putCoordGetIndices(
1412 double X1, double X2, double X3, int &Index1, int &Index2, int &Index3) {
1413 double pnt[3] = {X1, X2, X3};
1414 int index[3] = {0, 0, 0};
1415
1416 if (GridType < 5) {
1417 for (int i = 0; i < 3; ++i) {
1418 if (EasyCoordinate[i]) {
1419 index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
1420 }
1421 }
1422 for (int i = 0; i < 3; ++i) {
1423 if (!EasyCoordinate[i]) {
1424 double stepSize = BasicDistance0[i];
1425 double offset = 0.0;
1426 for (int j = 0; j < 3; ++j) {
1427 stepSize += BasicDistance1[i][j] * index[j];
1428 offset += BasicDistance2[i][j] * index[j];
1429 }
1430 index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
1431 }
1432 }
1433 }
1434 if (GridType == 5 || GridType == 6) {
1435 double sinPhi;
1436 if (GridType == 6) {
1437 sinPhi = cos(pnt[1]);
1438 } else if (GridType == 5) {
1439 sinPhi = sin(pnt[1]);
1440 } else {
1441 abort();
1442 }
1443
1444 if (std::abs(sinPhi) < EPSILON) {
1445 sinPhi = EPSILON;
1446 cout << "ERROR DIVISION BY ZERO" << endl;
1447 }
1448 double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1449 stepSize = stepSize / (NumberOfPoints[0] - 1);
1450 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1451 index[0] = int((pnt[0] - startingPoint) / stepSize);
1452 index[1] = int((pnt[1] - ReferencePoint[1]) / BasicDistance0[1]);
1453 index[2] = int((pnt[2] - ReferencePoint[2]) / BasicDistance0[2]);
1454 }
1455
1456 Index1 = index[0];
1457 Index2 = index[1];
1458 Index3 = index[2];
1459
1460 return;
1461 }
1462
1463 void prepareMagneticFieldGrid::putIndicesGetXAndB(
1464 int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1465 SixDPoint GridPoint;
1466 GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1467 X1 = GridPoint.x1();
1468 X2 = GridPoint.x2();
1469 X3 = GridPoint.x3();
1470 Bx = GridPoint.bx();
1471 By = GridPoint.by();
1472 Bz = GridPoint.bz();
1473
1474 return;
1475 }
1476
1477 void prepareMagneticFieldGrid::putIndCalcXReturnB(
1478 int Index1, int Index2, int Index3, double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1479 int index[3] = {Index1, Index2, Index3};
1480 double pnt[3] = {0., 0., 0.};
1481
1482 if (GridType < 5) {
1483 for (int i = 0; i < 3; ++i) {
1484 if (EasyCoordinate[i]) {
1485 pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
1486 } else {
1487 double stepSize = BasicDistance0[i];
1488 double offset = 0.0;
1489 for (int j = 0; j < 3; ++j) {
1490 stepSize += BasicDistance1[i][j] * index[j];
1491 offset += BasicDistance2[i][j] * index[j];
1492 }
1493 pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
1494 }
1495 }
1496 }
1497 if (GridType == 5 || GridType == 6) {
1498 pnt[2] = ReferencePoint[2] + BasicDistance0[2] * index[2];
1499 pnt[1] = ReferencePoint[1] + BasicDistance0[1] * index[1];
1500
1501 double sinPhi;
1502 if (GridType == 6) {
1503 sinPhi = cos(pnt[1]);
1504 } else if (GridType == 5) {
1505 sinPhi = sin(pnt[1]);
1506 } else {
1507 cout << "unknown GridType!" << endl;
1508 abort();
1509 }
1510
1511 if (std::abs(sinPhi) < EPSILON) {
1512 sinPhi = EPSILON;
1513 cout << "ERROR DIVISION BY ZERO" << endl;
1514 }
1515 double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
1516 stepSize = stepSize / (NumberOfPoints[0] - 1);
1517 double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
1518 pnt[0] = startingPoint + stepSize * index[0];
1519 }
1520
1521 X1 = pnt[0];
1522 X2 = pnt[1];
1523 X3 = pnt[2];
1524
1525 SixDPoint GridPoint;
1526 GridPoint = GridData.operator[](lineNumber(Index1, Index2, Index3));
1527 Bx = GridPoint.bx();
1528 By = GridPoint.by();
1529 Bz = GridPoint.bz();
1530
1531 return;
1532 }
1533
1534 int prepareMagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3) {
1535 return Index1 * NumberOfPoints[1] * NumberOfPoints[2] + Index2 * NumberOfPoints[2] + Index3;
1536 }
1537
1538 bool prepareMagneticFieldGrid::IndexedDoubleVector::operator<(const IndexedDoubleVector &x) const {
1539 if (I3[0] < x.I3[0])
1540 return true;
1541 else if (I3[0] == x.I3[0]) {
1542 if (I3[1] < x.I3[1])
1543 return true;
1544 else if (I3[1] == x.I3[1]) {
1545 if (I3[2] < x.I3[2])
1546 return true;
1547 else
1548 return false;
1549 } else
1550 return false;
1551 } else
1552 return false;
1553 }
1554
1555 void prepareMagneticFieldGrid::IndexedDoubleVector::putI3(int index1, int index2, int index3) {
1556 I3[0] = index1;
1557 I3[1] = index2;
1558 I3[2] = index3;
1559 return;
1560 }
1561
1562 void prepareMagneticFieldGrid::IndexedDoubleVector::getI3(int &index1, int &index2, int &index3) {
1563 index1 = I3[0];
1564 index2 = I3[1];
1565 index3 = I3[2];
1566 return;
1567 }
1568
1569 void prepareMagneticFieldGrid::IndexedDoubleVector::putV6(
1570 double X1, double X2, double X3, double Bx, double By, double Bz) {
1571 V6[0] = X1;
1572 V6[1] = X2;
1573 V6[2] = X3;
1574 V6[3] = Bx;
1575 V6[4] = By;
1576 V6[5] = Bz;
1577 return;
1578 }
1579
1580 void prepareMagneticFieldGrid::IndexedDoubleVector::getV6(
1581 double &X1, double &X2, double &X3, double &Bx, double &By, double &Bz) {
1582 X1 = V6[0];
1583 X2 = V6[1];
1584 X3 = V6[2];
1585 Bx = V6[3];
1586 By = V6[4];
1587 Bz = V6[5];
1588 return;
1589 }
1590
1591 void prepareMagneticFieldGrid::SixDPoint::putP6(double X1, double X2, double X3, double Bx, double By, double Bz) {
1592 P6[0] = X1;
1593 P6[1] = X2;
1594 P6[2] = X3;
1595 P6[3] = Bx;
1596 P6[4] = By;
1597 P6[5] = Bz;
1598 return;
1599 }
1600
1601 double prepareMagneticFieldGrid::SixDPoint::x1() { return P6[0]; }
1602
1603 double prepareMagneticFieldGrid::SixDPoint::x2() { return P6[1]; }
1604
1605 double prepareMagneticFieldGrid::SixDPoint::x3() { return P6[2]; }
1606
1607 double prepareMagneticFieldGrid::SixDPoint::bx() { return P6[3]; }
1608
1609 double prepareMagneticFieldGrid::SixDPoint::by() { return P6[4]; }
1610
1611 double prepareMagneticFieldGrid::SixDPoint::bz() { return P6[5]; }