EgammaPhotonAlgos/src/ConversionBarrelEstimator.cc

0001 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0002 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0005 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0006 #include "TrackingTools/DetLayers/interface/rangesIntersect.h"
0007 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0008
0009 // zero value indicates incompatible ts - hit pair
0010 std::pair<bool, double> ConversionBarrelEstimator::estimate(const TrajectoryStateOnSurface& ts,
0011                                                             const TrackingRecHit& hit) const {
0012   std::pair<bool, double> result;
0013
0014   //std::cout << "  ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const TransientTrackingRecHit& hit) " << std::endl;
0015
0016   float tsPhi = ts.globalParameters().position().phi();
0017   GlobalPoint gp = hit.globalPosition();
0018   float rhPhi = gp.phi();
0019
0020   // allow a z fudge of 2 sigma
0021   float dz = 2. * sqrt(hit.localPositionError().yy());
0022   float zDiff = ts.globalParameters().position().z() - gp.z();
0023   float phiDiff = tsPhi - rhPhi;
0024   if (phiDiff > pi)
0025     phiDiff -= twopi;
0026   if (phiDiff < -pi)
0027     phiDiff += twopi;
0028
0029   // add the errors on the window and the point in quadrature
0030   float zrange = sqrt(theZRangeMax * theZRangeMax + dz * dz);
0031
0032   /*
0033   std::cout << "  BarrelEstimator ts local error " <<ts.localError().positionError()  << " hit local error " << hit.localPositionError() << std::endl;
0034   std::cout << "  BarrelEstimator:  RecHit at " << gp << " phi " << rhPhi << " eta " << gp.eta() <<  std::endl;
0035   std::cout << "  BarrelEstimator:  ts at " << ts.globalParameters().position() << " phi " <<ts.globalParameters().position().phi() << " eta " << ts.globalParameters().position().eta()<<  std::endl;
0036   std::cout << "                    zrange = +/-" << zrange << ", zDiff = " << zDiff << std::endl;
0037   std::cout << "                    thePhiRangeMin = " << thePhiRangeMin << ", thePhiRangeMax = " << thePhiRangeMax << ", phiDiff = " << phiDiff << std::endl;
0038   */
0039
0040   if (phiDiff < thePhiRangeMax && phiDiff > thePhiRangeMin && zDiff < zrange && zDiff > -zrange) {
0041     //    std::cout << "      estimator returns 1 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
0042     //    << thePhiRangeMax << " and zDiff " << zDiff << " < " << zrange << std::endl;
0043     // std::cout << " YES " << phiDiff << " " << zDiff << std::endl;
0044     // std::cout << "                  => RECHIT ACCEPTED " << std::endl;
0045
0046     result.first = true;
0047     result.second = phiDiff;
0048   } else {
0049     //     std::cout << "      estimator returns NOT ACCEPTED  with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
0050     //<< thePhiRangeMax << " and zDiff " << zDiff << " < " << theZRangeMax+dz << std::endl;
0051
0052     result.first = false;
0053     result.second = 0;
0054   }
0055
0056   return result;
0057 }
0058
0059 bool ConversionBarrelEstimator::estimate(const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const {
0060   typedef std::pair<float, float> Range;
0061   //  std::cout << "  ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) " << std::endl;
0062
0063   GlobalPoint trajPos(ts.globalParameters().position());
0064   Range trajZRange(trajPos.z() - 2. * theZRangeMax, trajPos.z() + 2. * theZRangeMax);
0065   Range trajPhiRange(trajPos.phi() + thePhiRangeMin, trajPos.phi() + thePhiRangeMax);
0066
0067   if (rangesIntersect(trajZRange, plane.zSpan()) &&
0068       rangesIntersect(trajPhiRange, plane.phiSpan(), [](auto x, auto y) { return Geom::phiLess(x, y); })) {
0069     //     std::cout << "   ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane)  IN RANGE " << std::endl;
0070     return true;
0071
0072   } else {
0073     //    std::cout << "   ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) NOT IN RANGE " << std::endl;
0074     return false;
0075   }
0076 }
0077
0078 MeasurementEstimator::Local2DVector ConversionBarrelEstimator::maximalLocalDisplacement(
0079     const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const {
0080   /*
0081   if ( ts.hasError() ) {
0082     LocalError le = ts.localError().positionError();
0083     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;
0084     return Local2DVector( sqrt(le.xx())*nSigmaCut(), sqrt(le.yy())*nSigmaCut());
0085   }
0086
0087   else return Local2DVector(9999,9999);
0088   */
0089   return Local2DVector(9999, 9999);
0090 }