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" 

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

}

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;
    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));

}