File indexing completed on 2024-04-06 12:22:31
0001
0002
0003
0004
0005
0006
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
0030
0031 void MagGeometryExerciser::testFindVolume(int ntry) {
0032 cout << endl << "-----------------------------------------------------" << endl << " findVolume(random) test" << endl;
0033
0034
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
0060 bool MagGeometryExerciser::testFindVolume(const GlobalPoint& gp) {
0061 float tolerance = 0.;
0062
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
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
0085 void MagGeometryExerciser::testInside(int ntry, float tolerance) {
0086 cout << "-----------------------------------------------------" << endl << " inside(random) test" << endl;
0087
0088
0089
0090
0091 if (false) {
0092 cout << "Known points:" << endl;
0093 testInside(GlobalPoint(0., 0., 0.));
0094 testInside(GlobalPoint(331.358, -278.042, -648.788));
0095 testInside(GlobalPoint(-426.188, 75.1483, -533.075));
0096 testInside(GlobalPoint(-586.27, 157.094, 150.702));
0097 testInside(GlobalPoint(-465.562, -465.616, -222.224));
0098 testInside(GlobalPoint(637.078, 170.737, -222.632));
0099 testInside(GlobalPoint(162.971, -608.294, -306.791));
0100 testInside(GlobalPoint(-633.925, -169.839, -390.589));
0101 testInside(GlobalPoint(170.358, 635.857, -535.735));
0102 testInside(GlobalPoint(166.836, 622.58, -513.872));
0103 testInside(GlobalPoint(-12.7086, -3.99366, -1313.01));
0104 testInside(GlobalPoint(106.178, -538.867, -69.2348));
0105 testInside(GlobalPoint(19.1496, 522.831, -1333));
0106 testInside(GlobalPoint(355.516, -479.729, -1333.01));
0107 testInside(GlobalPoint(157.96, -389.223, -1333.01));
0108 testInside(GlobalPoint(-464.338, 464.234, -473.616));
0109 testInside(GlobalPoint(17.474, -669.279, -1333.01));
0110 testInside(GlobalPoint(-453.683, -453.601, -204.895));
0111 testInside(GlobalPoint(165.092, -615.998, -494.625));
0112 testInside(GlobalPoint(-586.27, 157.094, -617.882));
0113 testInside(GlobalPoint(-142.226, 390.763, -553.809));
0114 testInside(GlobalPoint(-141.701, 389.32, -142.088));
0115 }
0116
0117
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
0126
0127
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
0139
0140 bool MagGeometryExerciser::testInside(const GlobalPoint& gp, float tolerance) {
0141 bool reportSuccess = false;
0142
0143 vector<MagVolume6Faces const*>& vols = volumes;
0144
0145
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
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;
0197 }
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208