test/stubs/VolumeGridTester.cc

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) {  // move x by ulp times the float numerical precision
0014     return x + std::numeric_limits<float>::epsilon() * std::fabs(x) * ulp;
0015   }
0016 }  // namespace
0017
0018 bool VolumeGridTester::testInside() const {
0019   //   static string lastName("firstcall");
0020   //   if (lastName == volume_->name) return true; // skip multiple calls
0021   //   else lastName = volume_->name;
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   //   cout << "The volume position is " << volume_->position() << endl;
0033   //   cout << "Is the volume position inside the volume? " << volume_->inside(volume_->position(), tolerance) << endl;
0034
0035   Dimensions sizes = grid->dimensions();
0036   //   cout << "Grid has " << 3 << " dimensions "
0037   //        << " number of nodes is " << sizes.w << " " << sizes.h << " " << sizes.d << endl;
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         // Check that grid point is inside its own volume
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         // Check that points within numerical accuracy to grid boundary points are discoverable
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 }