File indexing completed on 2024-04-06 12:22:31
0001 #include "VolumeGridTester.h"
0002 #include "MagneticField/Interpolation/interface/MFGrid.h"
0003 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0004 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0005 #include <iostream>
0006 #include <iomanip>
0007 #include <string>
0008 #include <map>
0009
0010 using namespace std;
0011
0012 namespace {
0013 float eps(float x, int ulp) {
0014 return x + std::numeric_limits<float>::epsilon() * std::fabs(x) * ulp;
0015 }
0016 }
0017
0018 bool VolumeGridTester::testInside() const {
0019
0020
0021
0022
0023 const MFGrid* grid = dynamic_cast<const MFGrid*>(magProvider_);
0024 if (grid == 0) {
0025 cout << "VolumeGridTester: magProvider is not a MFGrid3D, cannot test it..." << endl << "expected ";
0026 return false;
0027 }
0028
0029 bool result = true;
0030 const double tolerance = 0.03;
0031
0032
0033
0034
0035 Dimensions sizes = grid->dimensions();
0036
0037
0038
0039 size_t dumpCount = 0;
0040 for (int j = 0; j < sizes.h; j++) {
0041 for (int k = 0; k < sizes.d; k++) {
0042 for (int i = 0; i < sizes.w; i++) {
0043 MFGrid::LocalPoint lp = grid->nodePosition(i, j, k);
0044
0045 if (!volume_->inside(lp, tolerance)) {
0046 result = false;
0047 if (++dumpCount < 2)
0048 dumpProblem(lp, tolerance);
0049 else
0050 return result;
0051 }
0052
0053
0054 if (field_ != nullptr &&
0055 ((j == 0 || j == sizes.h - 1) || (k == 0 || k == sizes.d - 1) || (i == 0 || i == sizes.w - 1))) {
0056 GlobalPoint gp = volume_->toGlobal(lp);
0057
0058 for (int nulp = -5; nulp <= 5; ++nulp) {
0059 result &= testFind(GlobalPoint(eps(gp.x(), nulp), gp.y(), gp.z()));
0060 result &= testFind(GlobalPoint(gp.x(), eps(gp.y(), nulp), gp.z()));
0061 result &= testFind(GlobalPoint(gp.x(), gp.y(), eps(gp.z(), nulp)));
0062
0063 if (result == false)
0064 return result;
0065 }
0066 }
0067 }
0068 }
0069 }
0070 return result;
0071 }
0072
0073 void VolumeGridTester::dumpProblem(const MFGrid::LocalPoint& lp, double tolerance) const {
0074 MFGrid::GlobalPoint gp(volume_->toGlobal(lp));
0075 cout << "ERROR: VolumeGridTester: Grid point " << lp << " (local) " << gp << " (global) " << gp.perp() << " "
0076 << gp.phi() << " (R,phi global) is not in its own volume!" << endl;
0077
0078 const vector<VolumeSide>& faces = volume_->faces();
0079 for (vector<VolumeSide>::const_iterator v = faces.begin(); v != faces.end(); v++) {
0080 cout << " Volume face has position " << v->surface().position() << " side " << (int)v->surfaceSide() << " rotation "
0081 << endl
0082 << v->surface().rotation() << endl;
0083
0084 Surface::Side side = v->surface().side(gp, tolerance);
0085 if (side != v->surfaceSide() && side != SurfaceOrientation::onSurface) {
0086 cout << " Wrong side: " << (int)side << " local position in surface frame " << v->surface().toLocal(gp) << endl;
0087 } else
0088 cout << " Correct side: " << (int)side << endl;
0089 }
0090 }
0091
0092 bool VolumeGridTester::testFind(GlobalPoint gp) const {
0093 if (field_->isDefined(gp)) {
0094 MagVolume const* vol = field_->findVolume(gp);
0095 if (vol == nullptr) {
0096 cout << "ERROR: VolumeGridTester: No volume found! Global point: " << setprecision(8) << gp
0097 << " , at R= " << gp.perp() << ", Z= " << gp.z() << endl;
0098 return false;
0099 }
0100 }
0101 return true;
0102 }