00001 #include "CLHEP/Units/PhysicalConstants.h" 00002 #include "FWCore/MessageLogger/interface/MessageLogger.h" 00003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h" 00004 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" 00005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" 00006 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h" 00007 #include "DataFormats/GeometrySurface/interface/BoundPlane.h" 00008 #include "TrackingTools/DetLayers/interface/PhiLess.h" 00009 #include "TrackingTools/DetLayers/interface/rangesIntersect.h" 00010 #include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h" 00011 00012 00013 00014 00015 // zero value indicates incompatible ts - hit pair 00016 std::pair<bool,double> ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts, 00017 const TransientTrackingRecHit& hit) const { 00018 LogDebug("ConversionForwardEstimator") << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts ...) " << "\n"; 00019 // std::cout << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts ...) " << "\n"; 00020 00021 std::pair<bool,double> result; 00022 00023 float tsPhi = ts.globalParameters().position().phi(); 00024 GlobalPoint gp = hit.globalPosition(); 00025 float rhPhi = gp.phi(); 00026 float rhR = gp.perp(); 00027 00028 // allow an r fudge of 1.5 * times the sigma 00029 // nodt used float dr = 1.5 * hit.localPositionError().yy(); 00030 //std::cout << " err " << hit.globalPositionError().phierr(gp) 00031 // << " " << hit.globalPositionError().rerr(gp) << std::endl; 00032 00033 // not used float zLayer = ts.globalParameters().position().z(); 00034 float rLayer = ts.globalParameters().position().perp(); 00035 00036 float newdr = sqrt(pow(dr_,2)+4.*hit.localPositionError().yy()); 00037 float rMin = rLayer - newdr; 00038 float rMax = rLayer + newdr; 00039 float phiDiff = tsPhi - rhPhi; 00040 if (phiDiff > pi) phiDiff -= twopi; 00041 if (phiDiff < -pi) phiDiff += twopi; 00042 00043 //std::cout << " ConversionForwardEstimator: RecHit at " << gp << "\n"; 00044 //std::cout << " rMin = " << rMin << ", rMax = " << rMax << ", rHit = " << rhR << "\n"; 00045 //std::cout << " thePhiRangeMin = " << thePhiRangeMin << ", thePhiRangeMax = " << thePhiRangeMax << ", phiDiff = " << phiDiff << "\n"; 00046 00047 00048 if ( phiDiff < thePhiRangeMax && phiDiff > thePhiRangeMin && 00049 rhR < rMax && rhR > rMin) { 00050 00051 00052 // std::cout << " estimator returns 1 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < " 00053 // << thePhiRangeMax << " and rhR " << rMin << " < " << rhR << " < " << rMax << "\n"; 00054 //std::cout << " YES " << phiDiff << " " <<rLayer-rhR << "\n"; 00055 //std::cout << " => RECHIT ACCEPTED " << "\n"; 00056 00057 result.first= true; 00058 result.second=phiDiff; 00059 } else { 00060 /* 00061 cout << " estimator returns 0 with phiDiff " << thePhiRangeMin << " < " << phiDiff << " < " 00062 << thePhiRangeMax << " and rhR " << rMin << " < " << rhR << " < " << rMax << endl; 00063 */ 00064 result.first= false; 00065 result.second=0; 00066 00067 } 00068 00069 return result; 00070 00071 } 00072 00073 bool ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts, 00074 const BoundPlane& plane) const { 00075 00076 // std::cout << "ConversionForwardEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) always TRUE " << "\n"; 00077 // this method should return one if a detector ring is close enough 00078 // to the hit, zero otherwise. 00079 // Now time is wasted looking for hits in the rings which are anyhow 00080 // too far from the prediction 00081 return true ; 00082 00083 } 00084 00085 00086 00087 MeasurementEstimator::Local2DVector 00088 ConversionForwardEstimator::maximalLocalDisplacement( const TrajectoryStateOnSurface& ts, 00089 const BoundPlane& plane) const 00090 { 00091 00092 /* 00093 if ( ts.hasError() ) { 00094 LocalError le = ts.localError().positionError(); 00095 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; 00096 return Local2DVector( sqrt(le.xx())*nSigmaCut(), sqrt(le.yy())*nSigmaCut()); 00097 00098 } 00099 else return Local2DVector(99999,99999); 00100 */ 00101 00102 return Local2DVector(99999,99999); 00103 00104 } 00105 00106