ConversionBarrelEstimator.cc

Go to the documentation of this file.

#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" 
#include "TrackingTools/DetLayers/interface/rangesIntersect.h"
#include "DataFormats/GeometryVector/interface/VectorUtil.h"


  // zero value indicates incompatible ts - hit pair
std::pair<bool,double> ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, 
                                const TrackingRecHit& hit) const {
  std::pair<bool,double> result;
  
  //std::cout << "  ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const TransientTrackingRecHit& hit) " << std::endl;
 
  float tsPhi = ts.globalParameters().position().phi();
  GlobalPoint gp = hit.globalPosition();
  float rhPhi = gp.phi();
  
  // allow a z fudge of 2 sigma
  float dz = 2. * sqrt(hit.localPositionError().yy()) ;
  float zDiff = ts.globalParameters().position().z() - gp.z(); 
  float phiDiff = tsPhi - rhPhi;
  if (phiDiff > pi) phiDiff -= twopi;
  if (phiDiff < -pi) phiDiff += twopi; 
  
  // add the errors on the window and the point in quadrature
  float zrange = sqrt(theZRangeMax*theZRangeMax + dz*dz);
  
  /*
  std::cout << "  BarrelEstimator ts local error " <<ts.localError().positionError()  << " hit local error " << hit.localPositionError() << std::endl; 
  std::cout << "  BarrelEstimator:  RecHit at " << gp << " phi " << rhPhi << " eta " << gp.eta() <<  std::endl;
  std::cout << "  BarrelEstimator:  ts at " << ts.globalParameters().position() << " phi " <<ts.globalParameters().position().phi() << " eta " << ts.globalParameters().position().eta()<<  std::endl;
  std::cout << "                    zrange = +/-" << zrange << ", zDiff = " << zDiff << std::endl;
  std::cout << "                    thePhiRangeMin = " << thePhiRangeMin << ", thePhiRangeMax = " << thePhiRangeMax << ", phiDiff = " << phiDiff << std::endl;
  */

  
  
  if ( phiDiff < thePhiRangeMax && phiDiff > thePhiRangeMin && 
       zDiff < zrange && zDiff > -zrange) {
    
    //    std::cout << "      estimator returns 1 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
    //    << thePhiRangeMax << " and zDiff " << zDiff << " < " << zrange << std::endl;
    // std::cout << " YES " << phiDiff << " " << zDiff << std::endl;
    // std::cout << "                  => RECHIT ACCEPTED " << std::endl;
    
    result.first=true;
    result.second=phiDiff;
  } else {
    
    //     std::cout << "      estimator returns NOT ACCEPTED  with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
    //<< thePhiRangeMax << " and zDiff " << zDiff << " < " << theZRangeMax+dz << std::endl;
    
    result.first=false;
    result.second=0;
    
  }
  
  return result;
  
}



bool ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, 
                               const BoundPlane& plane) const {
  
  typedef     std::pair<float,float>   Range;
  //  std::cout << "  ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) " << std::endl;

  GlobalPoint trajPos(ts.globalParameters().position());
  Range trajZRange(trajPos.z() - 2.*theZRangeMax, trajPos.z() + 2.*theZRangeMax);
  Range trajPhiRange(trajPos.phi() + thePhiRangeMin, trajPos.phi() + thePhiRangeMax);


   if(rangesIntersect(trajZRange, plane.zSpan()) &&
      rangesIntersect(trajPhiRange, plane.phiSpan(), [](auto x, auto y){return Geom::phiLess(x, y);})) {
     //     std::cout << "   ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane)  IN RANGE " << std::endl;  
    return true;   



  } else {

    //    std::cout << "   ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) NOT IN RANGE " << std::endl;  
    return false;

  }


}


MeasurementEstimator::Local2DVector
ConversionBarrelEstimator::maximalLocalDisplacement( const TrajectoryStateOnSurface& ts,
                                                        const BoundPlane& plane) const
{
  

  
  /* 
  if ( ts.hasError() ) {
    LocalError le = ts.localError().positionError();
    std::cout << "  ConversionBarrelEstimator::maximalLocalDisplacent local error " << sqrt(le.xx()) << " " << sqrt(le.yy()) << " nSigma " << nSigmaCut() << " sqrt(le.xx())*nSigmaCut() " << sqrt(le.xx())*nSigmaCut()  << "  sqrt(le.yy())*nSigmaCut() " <<  sqrt(le.yy())*nSigmaCut() << std::endl;
    return Local2DVector( sqrt(le.xx())*nSigmaCut(), sqrt(le.yy())*nSigmaCut());
  }

  else return Local2DVector(9999,9999);
  */
return Local2DVector(9999,9999);

}