test/stubs/MagGeometryExerciser.cc

0001 //#include "Utilities/Configuration/interface/Architecture.h"
0002
0003 /*
0004  *  See header file for a description of this class.
0005  *
0006  *  \author N. Amapane - INFN Torino
0007  */
0008
0009 #include "MagneticField/GeomBuilder/test/stubs/MagGeometryExerciser.h"
0010 #include "MagneticField/VolumeBasedEngine/interface/MagGeometry.h"
0011 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0012 #include "GlobalPointProvider.h"
0013
0014 #include <algorithm>
0015
0016 using namespace std;
0017
0018 MagGeometryExerciser::MagGeometryExerciser(const MagGeometry* g) : theGeometry(g) {
0019   const vector<MagVolume6Faces const*>& theBVolumes = theGeometry->barrelVolumes();
0020   const vector<MagVolume6Faces const*>& theEVolumes = theGeometry->endcapVolumes();
0021
0022   volumes = theBVolumes;
0023   volumes.insert(volumes.end(), theEVolumes.begin(), theEVolumes.end());
0024 }
0025
0026 MagGeometryExerciser::~MagGeometryExerciser() {}
0027
0028 //----------------------------------------------------------------------
0029 // Check if findVolume succeeds for random points.
0030 // Note: findVolumeTolerance in pset to change the tolerance.
0031 void MagGeometryExerciser::testFindVolume(int ntry) {
0032   cout << endl << "-----------------------------------------------------" << endl << " findVolume(random) test" << endl;
0033
0034   // Test  known overlaps/gaps
0035   if (true) {
0036     cout << "Known points:" << endl;
0037     testFindVolume(GlobalPoint(0, 0, 0));
0038   }
0039
0040   float maxZ = 2000.;
0041   if (theGeometry->geometryVersion >= 160812)
0042     maxZ = 2400.;
0043
0044   GlobalPointProvider p(0., 900., -Geom::pi(), Geom::pi(), -maxZ, maxZ);
0045
0046   cout << "Random points: (|Z| < " << maxZ << ")" << endl;
0047   int success = 0;
0048   for (int i = 0; i < ntry; ++i) {
0049     if (testFindVolume(p.getPoint())) {
0050       ++success;
0051     }
0052   }
0053
0054   cout << " Tested " << ntry << " Failures: " << ntry - success << endl
0055        << "-----------------------------------------------------" << endl;
0056 }
0057
0058 //----------------------------------------------------------------------
0059 // Check if findVolume succeeds for the given point.
0060 bool MagGeometryExerciser::testFindVolume(const GlobalPoint& gp) {
0061   float tolerance = 0.;
0062   //  float tolerance = 0.03;  // Note: findVolume should handle tolerance himself.
0063   MagVolume6Faces const* vol = (MagVolume6Faces const*)theGeometry->findVolume(gp, tolerance);
0064   bool ok = (vol != 0);
0065
0066   if (vol == 0) {
0067     cout << "Test ERROR: No volume found! Global point: " << gp << " , GP Z = " << gp.z() << ", GP perp = " << gp.perp()
0068          << " isBarrel: " << theGeometry->inBarrel(gp) << endl;
0069
0070     // Try with a linear search
0071     vol = (MagVolume6Faces const*)theGeometry->findVolume1(gp, tolerance);
0072     cout << "Was in volume: ";
0073     if (vol != 0)
0074       cout << vol->volumeNo << ":" << int(vol->copyno);
0075     else
0076       cout << "-1";
0077     cout << endl;
0078   }
0079
0080   return ok;
0081 }
0082
0083 //----------------------------------------------------------------------
0084 // Check that a set of points is inside() one and only one volume.
0085 void MagGeometryExerciser::testInside(int ntry, float tolerance) {
0086   cout << "-----------------------------------------------------" << endl << " inside(random) test" << endl;
0087
0088   // Test random points: they should be found inside() one and only one volume
0089
0090   // Test some known overlaps/gaps
0091   if (false) {  //FIXME
0092     cout << "Known points:" << endl;
0093     testInside(GlobalPoint(0., 0., 0.));
0094     testInside(GlobalPoint(331.358, -278.042, -648.788));   //V_ZN_81 V_ZN_82
0095     testInside(GlobalPoint(-426.188, 75.1483, -533.075));   //V_ZN_81 V_ZN_82
0096     testInside(GlobalPoint(-586.27, 157.094, 150.702));     //V_ZN_170 V_ZN_174
0097     testInside(GlobalPoint(-465.562, -465.616, -222.224));  //203 205
0098     testInside(GlobalPoint(637.078, 170.737, -222.632));    //203 205
0099     testInside(GlobalPoint(162.971, -608.294, -306.791));   //203 205
0100     testInside(GlobalPoint(-633.925, -169.839, -390.589));  //203 205
0101     testInside(GlobalPoint(170.358, 635.857, -535.735));    //207 209
0102     testInside(GlobalPoint(166.836, 622.58, -513.872));     //207 209
0103     testInside(GlobalPoint(-12.7086, -3.99366, -1313.01));  //53 56
0104     testInside(GlobalPoint(106.178, -538.867, -69.2348));   //147 148
0105     testInside(GlobalPoint(19.1496, 522.831, -1333));       //50 64
0106     testInside(GlobalPoint(355.516, -479.729, -1333.01));   //50 64
0107     testInside(GlobalPoint(157.96, -389.223, -1333.01));    //50 64
0108     testInside(GlobalPoint(-464.338, 464.234, -473.616));   // 207 209
0109     testInside(GlobalPoint(17.474, -669.279, -1333.01));    // 50 64
0110     testInside(GlobalPoint(-453.683, -453.601, -204.895));  // 203 205
0111     testInside(GlobalPoint(165.092, -615.998, -494.625));   // 207 209
0112     testInside(GlobalPoint(-586.27, 157.094, -617.882));    // 177 177
0113     testInside(GlobalPoint(-142.226, 390.763, -553.809));   //81 79
0114     testInside(GlobalPoint(-141.701, 389.32, -142.088));    //81 79
0115   }
0116
0117   // Full CMS
0118
0119   float maxZ = 1999.9;
0120   if (theGeometry->geometryVersion >= 160812)
0121     maxZ = 2399.9;
0122
0123   GlobalPointProvider p(0, 900, -Geom::pi(), Geom::pi(), -maxZ, maxZ);
0124
0125   // Zoom of one sector
0126   //  GlobalPointProvider p(350.,900.,-0.27,0.27,-1999.9,1999.9);
0127   //  GlobalPointProvider p(0.,900.,-0.27,0.27,-1999.9,1999.9);
0128
0129   cout << "Random points:" << volumes.size() << " volumes" << endl;
0130   for (int i = 0; i < ntry; ++i) {
0131     if (i % 1000 == 0)
0132       cout << "test # " << i << endl;
0133     testInside(p.getPoint(), tolerance);
0134   }
0135 }
0136
0137 //----------------------------------------------------------------------
0138 // Check that the given point is inside() one and only one volume.
0139 // This is not always the case due to thin gaps and tolerance...
0140 bool MagGeometryExerciser::testInside(const GlobalPoint& gp, float tolerance) {
0141   bool reportSuccess = false;
0142
0143   vector<MagVolume6Faces const*>& vols = volumes;
0144   // or use only barrel volumes:
0145   // const vector<MagVolume6Faces const*>& vols = theGeometry->barrelVolumes();
0146
0147   MagVolume6Faces const* found = 0;
0148   for (vector<MagVolume6Faces const*>::const_iterator v = vols.begin(); v != vols.end(); ++v) {
0149     if ((*v) == 0) {
0150       cout << endl << "ERROR: no magvlolume" << endl;
0151       continue;
0152     }
0153     if ((*v)->inside(gp, tolerance)) {
0154       if (reportSuccess)
0155         cout << gp << " is inside vol: " << (*v)->volumeNo;
0156       if (found != 0) {
0157         cout << " ***ERROR: for " << gp << " found " << (*v)->volumeNo << " volume already found: " << found->volumeNo
0158              << endl;
0159       }
0160       found = (*v);
0161     }
0162   }
0163
0164   if (found == 0) {
0165     MagVolume6Faces const* foundP = 0;
0166     MagVolume6Faces const* foundN = 0;
0167     // Look for the closest neighbouring volumes
0168     const float phi = gp.phi();
0169     GlobalPoint gpP, gpN;
0170     int ntry = 0;
0171     while ((foundP == 0 || foundP == 0) && ntry < 60) {
0172       ++ntry;
0173       for (vector<MagVolume6Faces const*>::const_iterator v = vols.begin(); v != vols.end(); ++v) {
0174         if (foundP == 0) {
0175           float phiP = phi + ntry * 0.008727;
0176           GlobalPoint gpP(GlobalPoint::Cylindrical(gp.perp(), phiP, gp.z()));
0177           if ((*v)->inside(gpP)) {
0178             foundP = (*v);
0179           }
0180         }
0181         if (foundN == 0) {
0182           float phiN = phi - ntry * 0.008727;
0183           GlobalPoint gpN(GlobalPoint::Cylindrical(gp.perp(), phiN, gp.z()));
0184           if ((*v)->inside(gpN)) {
0185             foundN = (*v);
0186           }
0187         }
0188       }
0189     }
0190
0191     cout << gp << " ***ERROR no volume found! : closests: " << ((foundP == 0) ? -1 : foundP->volumeNo)
0192          << " at dphi: " << gpP.phi() - phi << " " << ((foundN == 0) ? -1 : foundN->volumeNo)
0193          << " at dphi: " << gpN.phi() - phi << endl;
0194   }
0195
0196   return true;  //FIXME
0197 }
0198
0199 //----------------------------------------------------------------------
0200
0201 //   if (test82) {
0202 //      minZ = -660;
0203 //      maxZ = 660.;
0204 //      minR = 411.5; //V81
0205 //      maxR = 447.; //V83
0206 //      minPhi= 104./180.*Geom::pi();
0207 //      maxPhi= 113./180.*Geom::pi();
0208 //   }