VolumeBasedEngine/src/MagGeometry.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino
0005  */
0006
0007 #include "MagneticField/VolumeBasedEngine/interface/MagGeometry.h"
0008 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
0009 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0010 #include "MagneticField/Layers/interface/MagBLayer.h"
0011 #include "MagneticField/Layers/interface/MagESector.h"
0012
0013 #include "Utilities/BinningTools/interface/PeriodicBinFinderInPhi.h"
0014
0015 #include "FWCore/Utilities/interface/isFinite.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include <iostream>
0019
0020 using namespace std;
0021 using namespace edm;
0022
0023 namespace {
0024   // A thread-local cache is accepted because MagneticField is used by several other ESProducts
0025   // via member variable, MagneticField and the other ESProducts are widely used, and migrating
0026   // all the uses of all those was deemed to have very high cost.
0027   std::atomic<int> instanceCounter(0);
0028   thread_local int localInstance = 0;
0029   thread_local MagVolume const* lastVolume = nullptr;
0030 }  // namespace
0031
0032 MagGeometry::MagGeometry(int geomVersion,
0033                          const std::vector<MagBLayer*>& tbl,
0034                          const std::vector<MagESector*>& tes,
0035                          const std::vector<MagVolume6Faces*>& tbv,
0036                          const std::vector<MagVolume6Faces*>& tev)
0037     : MagGeometry(geomVersion,
0038                   reinterpret_cast<std::vector<MagBLayer const*> const&>(tbl),
0039                   reinterpret_cast<std::vector<MagESector const*> const&>(tes),
0040                   reinterpret_cast<std::vector<MagVolume6Faces const*> const&>(tbv),
0041                   reinterpret_cast<std::vector<MagVolume6Faces const*> const&>(tev)) {}
0042
0043 MagGeometry::MagGeometry(int geomVersion,
0044                          const std::vector<MagBLayer const*>& tbl,
0045                          const std::vector<MagESector const*>& tes,
0046                          const std::vector<MagVolume6Faces const*>& tbv,
0047                          const std::vector<MagVolume6Faces const*>& tev)
0048     : me_(++instanceCounter),
0049       theBLayers(tbl),
0050       theESectors(tes),
0051       theBVolumes(tbv),
0052       theEVolumes(tev),
0053       cacheLastVolume(true),
0054       geometryVersion(geomVersion) {
0055   vector<double> rBorders;
0056
0057   for (vector<MagBLayer const*>::const_iterator ilay = theBLayers.begin(); ilay != theBLayers.end(); ++ilay) {
0058     LogTrace("MagGeoBuilder") << "  Barrel layer at " << (*ilay)->minR() << endl;
0059     //FIXME assume layers are already sorted in minR
0060     rBorders.push_back((*ilay)->minR() * (*ilay)->minR());
0061   }
0062
0063   theBarrelBinFinder = new MagBinFinders::GeneralBinFinderInR<double>(rBorders);
0064
0065 #ifdef EDM_ML_DEBUG
0066   for (vector<MagESector const*>::const_iterator isec = theESectors.begin(); isec != theESectors.end(); ++isec) {
0067     LogTrace("MagGeoBuilder") << "  Endcap sector at " << (*isec)->minPhi() << endl;
0068   }
0069 #endif
0070
0071   //FIXME assume sectors are already sorted in phi
0072   //FIXME: PeriodicBinFinderInPhi gets *center* of first bin
0073   int nEBins = theESectors.size();
0074   if (nEBins > 0)
0075     theEndcapBinFinder = new PeriodicBinFinderInPhi<float>(theESectors.front()->minPhi() + Geom::pi() / nEBins, nEBins);
0076
0077   // Compute barrel dimensions based on geometry version
0078   // FIXME: it would be nice to derive these from the actual geometry in the builder, possibly adding some specification to the geometry.
0079   switch (geomVersion >= 120812 ? 0 : (geomVersion >= 90812 ? 1 : 2)) {
0080     case 0:  // since 120812
0081       theBarrelRsq1 = 172.400 * 172.400;
0082       theBarrelRsq2 = 308.735 * 308.735;
0083       theBarrelZ0 = 350.000;
0084       theBarrelZ1 = 633.290;
0085       theBarrelZ2 = 662.010;
0086       break;
0087     case 1:  // version 90812 (no longer in use)
0088       theBarrelRsq1 = 172.400 * 172.400;
0089       theBarrelRsq2 = 308.755 * 308.755;
0090       theBarrelZ0 = 350.000;
0091       theBarrelZ1 = 633.890;
0092       theBarrelZ2 = 662.010;
0093       break;
0094     case 2:  // versions 71212, 90322
0095       theBarrelRsq1 = 172.400 * 172.400;
0096       theBarrelRsq2 = 308.755 * 308.755;
0097       theBarrelZ0 = 350.000;
0098       theBarrelZ1 = 633.290;
0099       theBarrelZ2 = 661.010;
0100       break;
0101   }
0102
0103   LogTrace("MagGeometry_cache") << "*** In MagGeometry ctor: me_=" << me_ << " instanceCounter=" << instanceCounter
0104                                 << endl;
0105 }
0106
0107 MagGeometry::~MagGeometry() {
0108   if (theBarrelBinFinder != nullptr)
0109     delete theBarrelBinFinder;
0110   if (theEndcapBinFinder != nullptr)
0111     delete theEndcapBinFinder;
0112
0113   for (vector<MagBLayer const*>::const_iterator ilay = theBLayers.begin(); ilay != theBLayers.end(); ++ilay) {
0114     delete (*ilay);
0115   }
0116
0117   for (vector<MagESector const*>::const_iterator ilay = theESectors.begin(); ilay != theESectors.end(); ++ilay) {
0118     delete (*ilay);
0119   }
0120 }
0121
0122 // Return field vector at the specified global point
0123 GlobalVector MagGeometry::fieldInTesla(const GlobalPoint& gp) const {
0124   MagVolume const* v = nullptr;
0125
0126   v = findVolume(gp);
0127   if (v != nullptr) {
0128     return v->fieldInTesla(gp);
0129   }
0130
0131   // Fall-back case: no volume found
0132
0133   if (edm::isNotFinite(gp.mag())) {
0134     LogWarning("MagneticField") << "Input value invalid (not a number): " << gp << endl;
0135
0136   } else {
0137     LogWarning("MagneticField") << "MagGeometry::fieldInTesla: failed to find volume for " << gp << endl;
0138   }
0139   return GlobalVector();
0140 }
0141
0142 // Linear search implementation (just for testing)
0143 MagVolume const* MagGeometry::findVolume1(const GlobalPoint& gp, double tolerance) const {
0144   MagVolume6Faces const* found = nullptr;
0145
0146   int errCnt = 0;
0147   if (inBarrel(gp)) {  // Barrel
0148     for (vector<MagVolume6Faces const*>::const_iterator v = theBVolumes.begin(); v != theBVolumes.end(); ++v) {
0149       if ((*v) == nullptr) {  //FIXME: remove this check
0150         LogError("MagGeometry") << endl << "***ERROR: MagGeometry::findVolume: MagVolume for barrel not set" << endl;
0151         ++errCnt;
0152         if (errCnt < 3)
0153           continue;
0154         else
0155           break;
0156       }
0157       if ((*v)->inside(gp, tolerance)) {
0158         found = (*v);
0159         break;
0160       }
0161     }
0162
0163   } else {  // Endcaps
0164     for (vector<MagVolume6Faces const*>::const_iterator v = theEVolumes.begin(); v != theEVolumes.end(); ++v) {
0165       if ((*v) == nullptr) {  //FIXME: remove this check
0166         LogError("MagGeometry") << endl << "***ERROR: MagGeometry::findVolume: MagVolume for endcap not set" << endl;
0167         ++errCnt;
0168         if (errCnt < 3)
0169           continue;
0170         else
0171           break;
0172       }
0173       if ((*v)->inside(gp, tolerance)) {
0174         found = (*v);
0175         break;
0176       }
0177     }
0178   }
0179
0180   return found;
0181 }
0182
0183 // Use hierarchical structure for fast lookup.
0184 MagVolume const* MagGeometry::findVolume(const GlobalPoint& gp, double tolerance) const {
0185   // Clear volume cache if this is a new instance
0186   if (me_ != localInstance) {
0187     LogTrace("MagGeometry_cache") << "*** In MagGeometry::findVolume resetting cache: me=" << me_
0188                                   << " localInstance=" << localInstance << endl;
0189     localInstance = me_;
0190     lastVolume = nullptr;
0191   }
0192
0193   if (lastVolume != nullptr && lastVolume->inside(gp)) {
0194     return lastVolume;
0195   }
0196
0197   MagVolume const* result = nullptr;
0198   if (inBarrel(gp)) {  // Barrel
0199     double aRsq = gp.perp2();
0200     int bin = theBarrelBinFinder->binIndex(aRsq);
0201
0202     // Search up to 3 layers inwards. This may happen for very thin layers.
0203     for (int bin1 = bin; bin1 >= max(0, bin - 3); --bin1) {
0204       LogTrace("MagGeometry") << "Trying layer at R " << theBLayers[bin1]->minR() << " " << sqrt(aRsq) << endl;
0205       result = theBLayers[bin1]->findVolume(gp, tolerance);
0206       LogTrace("MagGeometry") << "***In blayer " << bin1 - bin << " " << (result == nullptr ? " failed " : " OK ")
0207                               << endl;
0208       if (result != nullptr)
0209         break;
0210     }
0211
0212   } else {  // Endcaps
0213     Geom::Phi<float> phi = gp.phi();
0214     if (theEndcapBinFinder != nullptr && !theESectors.empty()) {
0215       int bin = theEndcapBinFinder->binIndex(phi);
0216       LogTrace("MagGeometry") << "Trying endcap sector at phi " << theESectors[bin]->minPhi() << " " << phi << endl;
0217       result = theESectors[bin]->findVolume(gp, tolerance);
0218       LogTrace("MagGeometry") << "***In guessed esector " << (result == nullptr ? " failed " : " OK ") << endl;
0219     } else
0220       edm::LogError("MagGeometry") << "Endcap empty";
0221   }
0222
0223   if (result == nullptr && tolerance < 0.0001) {
0224     // If search fails, retry with a 300 micron tolerance.
0225     // This is a hack for thin gaps on air-iron boundaries,
0226     // which will not be present anymore once surfaces are matched.
0227     LogTrace("MagGeometry") << "Increasing the tolerance to 0.03" << endl;
0228     result = findVolume(gp, 0.03);
0229   }
0230
0231   if (cacheLastVolume)
0232     lastVolume = result;
0233
0234   return result;
0235 }
0236
0237 bool MagGeometry::inBarrel(const GlobalPoint& gp) const {
0238   double aZ = fabs(gp.z());
0239   double aRsq = gp.perp2();
0240
0241   return ((aZ < theBarrelZ0) || (aZ < theBarrelZ1 && aRsq > theBarrelRsq1) ||
0242           (aZ < theBarrelZ2 && aRsq > theBarrelRsq2));
0243 }