EgammaPhotonAlgos/src/ConversionForwardEstimator.cc

0001 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
0004 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0006 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0007
0008 // zero value indicates incompatible ts - hit pair
0009 std::pair<bool, double> ConversionForwardEstimator::estimate(const TrajectoryStateOnSurface& ts,
0010                                                              const TrackingRecHit& hit) const {
0011   LogDebug("ConversionForwardEstimator")
0012       << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts ...) "
0013       << "\n";
0014   //  std::cout  << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts ...) " << "\n";
0015
0016   std::pair<bool, double> result;
0017
0018   float tsPhi = ts.globalParameters().position().phi();
0019   GlobalPoint gp = hit.globalPosition();
0020   float rhPhi = gp.phi();
0021   float rhR = gp.perp();
0022
0023   // allow an r fudge of 1.5 * times the sigma
0024   // nodt used float dr = 1.5 * hit.localPositionError().yy();
0025   //std::cout << " err " << hit.globalPositionError().phierr(gp)
0026   //    << " "     << hit.globalPositionError().rerr(gp) << std::endl;
0027
0028   // not used float zLayer = ts.globalParameters().position().z();
0029   float rLayer = ts.globalParameters().position().perp();
0030
0031   float newdr = sqrt(pow(dr_, 2) + 4. * hit.localPositionError().yy());
0032   float rMin = rLayer - newdr;
0033   float rMax = rLayer + newdr;
0034   float phiDiff = tsPhi - rhPhi;
0035   if (phiDiff > pi)
0036     phiDiff -= twopi;
0037   if (phiDiff < -pi)
0038     phiDiff += twopi;
0039
0040   //std::cout << " ConversionForwardEstimator: RecHit at " << gp << "\n";
0041   //std::cout << "                   rMin = " << rMin << ", rMax = " << rMax << ", rHit = " << rhR << "\n";
0042   //std::cout << "                   thePhiRangeMin = " << thePhiRangeMin << ", thePhiRangeMax = " << thePhiRangeMax << ", phiDiff = " << phiDiff << "\n";
0043
0044   if (phiDiff < thePhiRangeMax && phiDiff > thePhiRangeMin && rhR < rMax && rhR > rMin) {
0045     //    std::cout << "      estimator returns 1 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
0046     // << thePhiRangeMax << " and rhR " << rMin << " < " << rhR << " < " << rMax << "\n";
0047     //std::cout << " YES " << phiDiff << " " <<rLayer-rhR << "\n";
0048     //std::cout << "                  => RECHIT ACCEPTED " << "\n";
0049
0050     result.first = true;
0051     result.second = phiDiff;
0052   } else {
0053     /*
0054     cout << "      estimator returns 0 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < "
0055      << thePhiRangeMax << " and  rhR " << rMin << " < " << rhR << " < " << rMax << endl;
0056     */
0057     result.first = false;
0058     result.second = 0;
0059   }
0060
0061   return result;
0062 }
0063
0064 bool ConversionForwardEstimator::estimate(const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const {
0065   //  std::cout << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) always TRUE " << "\n";
0066   // this method should return one if a detector ring is close enough
0067   //     to the hit, zero otherwise.
0068   //     Now time is wasted looking for hits in the rings which are anyhow
0069   //     too far from the prediction
0070   return true;
0071 }
0072
0073 MeasurementEstimator::Local2DVector ConversionForwardEstimator::maximalLocalDisplacement(
0074     const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const {
0075   /*
0076   if ( ts.hasError() ) {
0077     LocalError le = ts.localError().positionError();
0078     std::cout << "  ConversionForwardEstimator::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;
0079     return Local2DVector( sqrt(le.xx())*nSigmaCut(), sqrt(le.yy())*nSigmaCut());
0080
0081   }
0082   else return Local2DVector(99999,99999);
0083   */
0084
0085   return Local2DVector(99999, 99999);
0086 }