GeomDetCompatibilityChecker.cc

Go to the documentation of this file.

// #define STAT_TSB

#ifdef STAT_TSB
#include <iostream>
#endif

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/thread_safety_macros.h"
#include "TrackingTools/DetLayers/interface/GeomDetCompatibilityChecker.h"
#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "FWCore/Utilities/interface/Likely.h"

namespace {

  struct Stat {
    struct Nop {
      Nop(int = 0) {}
      Nop& operator=(int) { return *this; }
      void operator++(int) {}
    };
#ifdef STAT_TSB
    using VAR = long long;
#else
    using VAR = Nop;
#endif
    VAR ntot = 0;
    VAR nf1 = 0;
    VAR nf2 = 0;

    VAR ns1 = 0;
    VAR ns2 = 0;
    VAR ns11 = 0;
    VAR ns21 = 0;

    VAR nth = 0;
    VAR nle = 0;

    //Geom checker     1337311   84696       634946      20369      259701        241266       18435         614128
    //    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
    //    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
    //    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
    //    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
    //    Geom checker 214,884,027  6,756,424 54,479,049  7,059,696  78,135,883 18443999 1124058 158,174,933 17153503
    //    Geom checker 453,155,905 14,054,554 79,733,432 14,837,002 163,414,609 0 0 324,629,999 0
#ifdef STAT_TSB
    ~Stat() {
      std::cout << "Geom checker " << ntot << ' ' << nf1 << ' ' << nf2 << ' ' << ns1 << ' ' << ns2 << ' ' << ns11 << ' '
                << ns21 << ' ' << nth << ' ' << nle << std::endl;
    }
#endif
  };

  CMS_THREAD_SAFE Stat stat;  // for production purpose it is thread safe

}  // namespace

std::pair<bool, TrajectoryStateOnSurface> GeomDetCompatibilityChecker::isCompatible(const GeomDet* theDet,
                                                                                    const TrajectoryStateOnSurface& tsos,
                                                                                    const Propagator& prop,
                                                                                    const MeasurementEstimator& est) {
  stat.ntot++;

  auto const sagCut = est.maxSagitta();
  auto const minTol2 = est.minTolerance2();

  // std::cout << "param " << sagCut << ' ' << std::sqrt(minTol2) << std::endl;

  /*
  auto err2 = tsos.curvilinearError().matrix()(3,3);
  auto largeErr = err2> 0.1*tolerance2;
  if (largeErr) stat.nle++; 
  */

  bool isIn = false;
  float sagitta = 99999999.0f;
  bool close = false;
  if LIKELY (sagCut > 0) {
    // linear approximation
    auto const& plane = theDet->specificSurface();
    StraightLinePlaneCrossing crossing(
        tsos.globalPosition().basicVector(), tsos.globalMomentum().basicVector(), prop.propagationDirection());
    auto path = crossing.pathLength(plane);
    isIn = path.first;
    if UNLIKELY (!path.first)
      stat.ns1++;
    else {
      auto gpos = GlobalPoint(crossing.position(path.second));
      auto tpath2 = (gpos - tsos.globalPosition()).perp2();
      // sagitta = d^2*c/2
      sagitta = 0.5f * std::abs(tpath2 * tsos.globalParameters().transverseCurvature());
      close = sagitta < sagCut;
      LogDebug("TkDetLayer") << "GeomDetCompatibilityChecker: sagitta " << sagitta << std::endl;
      if (close) {
        stat.nth++;
        auto pos = plane.toLocal(GlobalPoint(crossing.position(path.second)));
        // auto toll = LocalError(tolerance2,0,tolerance2);
        auto tollL2 = std::max(sagitta * sagitta, minTol2);
        auto toll = LocalError(tollL2, 0, tollL2);
        isIn = plane.bounds().inside(pos, toll);
        if (!isIn) {
          stat.ns2++;
          LogDebug("TkDetLayer") << "GeomDetCompatibilityChecker: not in " << pos << std::endl;
          return std::make_pair(false, TrajectoryStateOnSurface());
        }
      }
    }
  }

  // precise propagation
  TrajectoryStateOnSurface&& propSt = prop.propagate(tsos, theDet->specificSurface());
  if UNLIKELY (!propSt.isValid()) {
    stat.nf1++;
    return std::make_pair(false, std::move(propSt));
  }

  auto es = est.estimate(propSt, theDet->specificSurface());
  if (!es)
    stat.nf2++;
  if (close && (!isIn) && (!es))
    stat.ns11++;
  if (close && es && (!isIn)) {
    stat.ns21++;
  }  // std::cout << sagitta << std::endl;}
  return std::make_pair(es, std::move(propSt));
}