DetLayers/src/GeomDetCompatibilityChecker.cc

0001 // #define STAT_TSB
0002
0003 #ifdef STAT_TSB
0004 #include <iostream>
0005 #endif
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0009 #include "TrackingTools/DetLayers/interface/GeomDetCompatibilityChecker.h"
0010 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012 #include "FWCore/Utilities/interface/Likely.h"
0013
0014 namespace {
0015
0016   struct Stat {
0017     struct Nop {
0018       Nop(int = 0) {}
0019       Nop& operator=(int) { return *this; }
0020       void operator++(int) {}
0021     };
0022 #ifdef STAT_TSB
0023     using VAR = long long;
0024 #else
0025     using VAR = Nop;
0026 #endif
0027     VAR ntot = 0;
0028     VAR nf1 = 0;
0029     VAR nf2 = 0;
0030
0031     VAR ns1 = 0;
0032     VAR ns2 = 0;
0033     VAR ns11 = 0;
0034     VAR ns21 = 0;
0035
0036     VAR nth = 0;
0037     VAR nle = 0;
0038
0039     //Geom checker     1337311   84696       634946      20369      259701        241266       18435         614128
0040     //    Geom checker 124,567,704  3,862,821 36,055,605 3,799,127 29,825,229    4,573,316     320,063         75,420,840
0041     //    Geom checker 119,618,014  2,307,939 31,142,922 2,903,245 34,673,978    5,139,741     539,152         86,847,116 18,196,497
0042     //    Geom checker 125,554,439  3,431,348 31,900,589 3,706,531 37,272,039    5,160,257     1,670,236       90,573,031 19,505,412
0043     //    Geom checker 119,583,440  2,379,307 28,357,175 2,960,173 38,977,837    6,239,242       620,636       86,726,732  9,574,902
0044     //    Geom checker 214,884,027  6,756,424 54,479,049  7,059,696  78,135,883 18443999 1124058 158,174,933 17153503
0045     //    Geom checker 453,155,905 14,054,554 79,733,432 14,837,002 163,414,609 0 0 324,629,999 0
0046 #ifdef STAT_TSB
0047     ~Stat() {
0048       std::cout << "Geom checker " << ntot << ' ' << nf1 << ' ' << nf2 << ' ' << ns1 << ' ' << ns2 << ' ' << ns11 << ' '
0049                 << ns21 << ' ' << nth << ' ' << nle << std::endl;
0050     }
0051 #endif
0052   };
0053
0054   CMS_THREAD_SAFE Stat stat;  // for production purpose it is thread safe
0055
0056 }  // namespace
0057
0058 std::pair<bool, TrajectoryStateOnSurface> GeomDetCompatibilityChecker::isCompatible(const GeomDet* theDet,
0059                                                                                     const TrajectoryStateOnSurface& tsos,
0060                                                                                     const Propagator& prop,
0061                                                                                     const MeasurementEstimator& est) {
0062   stat.ntot++;
0063
0064   auto const sagCut = est.maxSagitta();
0065   auto const minTol2 = est.minTolerance2();
0066
0067   // std::cout << "param " << sagCut << ' ' << std::sqrt(minTol2) << std::endl;
0068
0069   /*
0070   auto err2 = tsos.curvilinearError().matrix()(3,3);
0071   auto largeErr = err2> 0.1*tolerance2;
0072   if (largeErr) stat.nle++;
0073   */
0074
0075   bool isIn = false;
0076   float sagitta = 99999999.0f;
0077   bool close = false;
0078   if LIKELY (sagCut > 0) {
0079     // linear approximation
0080     auto const& plane = theDet->specificSurface();
0081     StraightLinePlaneCrossing crossing(
0082         tsos.globalPosition().basicVector(), tsos.globalMomentum().basicVector(), prop.propagationDirection());
0083     auto path = crossing.pathLength(plane);
0084     isIn = path.first;
0085     if UNLIKELY (!path.first)
0086       stat.ns1++;
0087     else {
0088       auto gpos = GlobalPoint(crossing.position(path.second));
0089       auto tpath2 = (gpos - tsos.globalPosition()).perp2();
0090       // sagitta = d^2*c/2
0091       sagitta = 0.5f * std::abs(tpath2 * tsos.globalParameters().transverseCurvature());
0092       close = sagitta < sagCut;
0093       LogDebug("TkDetLayer") << "GeomDetCompatibilityChecker: sagitta " << sagitta << std::endl;
0094       if (close) {
0095         stat.nth++;
0096         auto pos = plane.toLocal(GlobalPoint(crossing.position(path.second)));
0097         // auto toll = LocalError(tolerance2,0,tolerance2);
0098         auto tollL2 = std::max(sagitta * sagitta, minTol2);
0099         auto toll = LocalError(tollL2, 0, tollL2);
0100         isIn = plane.bounds().inside(pos, toll);
0101         if (!isIn) {
0102           stat.ns2++;
0103           LogDebug("TkDetLayer") << "GeomDetCompatibilityChecker: not in " << pos << std::endl;
0104           return std::make_pair(false, TrajectoryStateOnSurface());
0105         }
0106       }
0107     }
0108   }
0109
0110   // precise propagation
0111   TrajectoryStateOnSurface&& propSt = prop.propagate(tsos, theDet->specificSurface());
0112   if UNLIKELY (!propSt.isValid()) {
0113     stat.nf1++;
0114     return std::make_pair(false, std::move(propSt));
0115   }
0116
0117   auto es = est.estimate(propSt, theDet->specificSurface());
0118   if (!es)
0119     stat.nf2++;
0120   if (close && (!isIn) && (!es))
0121     stat.ns11++;
0122   if (close && es && (!isIn)) {
0123     stat.ns21++;
0124   }  // std::cout << sagitta << std::endl;}
0125   return std::make_pair(es, std::move(propSt));
0126 }