diff --git a/DataFormats/GeometrySurface/src/BoundSpan.cc b/DataFormats/GeometrySurface/src/BoundSpan.cc index 4bc8f17d33b0c..38f53f5f08494 100644 --- a/DataFormats/GeometrySurface/src/BoundSpan.cc +++ b/DataFormats/GeometrySurface/src/BoundSpan.cc @@ -14,10 +14,10 @@ void BoundSpan::compute(Surface const & plane) { if (trapezoidalBounds) { std::array const & parameters = (*trapezoidalBounds).parameters(); - float hbotedge = parameters[0]; - float htopedge = parameters[1]; - float hapothem = parameters[3]; - float thickness = (*trapezoidalBounds).thickness(); + auto hbotedge = parameters[0]; + auto htopedge = parameters[1]; + auto hapothem = parameters[3]; + auto thickness = (*trapezoidalBounds).thickness(); corners[0] = plane.toGlobal( LocalPoint( -htopedge, hapothem, thickness/2)); corners[1] = plane.toGlobal( LocalPoint( htopedge, hapothem, thickness/2)); @@ -29,9 +29,9 @@ void BoundSpan::compute(Surface const & plane) { corners[7] = plane.toGlobal( LocalPoint( -hbotedge, -hapothem, -thickness/2)); }else if(rectangularBounds) { - float length = rectangularBounds->length(); - float width = rectangularBounds->width(); - float thickness = (*rectangularBounds).thickness(); + auto length = rectangularBounds->length(); + auto width = rectangularBounds->width(); + auto thickness = (*rectangularBounds).thickness(); corners[0] = plane.toGlobal( LocalPoint( -width/2, -length/2, thickness/2)); corners[1] = plane.toGlobal( LocalPoint( -width/2, +length/2, thickness/2)); @@ -47,16 +47,23 @@ void BoundSpan::compute(Surface const & plane) { float phimin = corners[0].barePhi(); float phimax = phimin; float zmin = corners[0].z(); float zmax = zmin; + float rmin = corners[0].perp2(); float rmax = rmin; for ( int i = 1; i < 8; i++ ) { - float cPhi = corners[i].barePhi(); + auto cPhi = corners[i].barePhi(); if ( Geom::phiLess( cPhi, phimin)) { phimin = cPhi; } if ( Geom::phiLess( phimax, cPhi)) { phimax = cPhi; } - float z = corners[i].z(); + auto z = corners[i].z(); if ( z < zmin) zmin = z; - if ( z > zmax) zmax = z; + if ( z > zmax) zmax = z; + auto r = corners[i].perp2(); + if ( r < rmin) rmin = r; + if ( r > rmax) rmax = r; + } m_zSpan.first = zmin; m_zSpan.second = zmax; + m_rSpan.first = std::sqrt(rmin); + m_rSpan.second = std::sqrt(rmax); m_phiSpan.first = phimin; m_phiSpan.second = phimax; } diff --git a/DataFormats/GeometryVector/interface/Pi.h b/DataFormats/GeometryVector/interface/Pi.h index f2f581aacf394..a7545df59f5cd 100644 --- a/DataFormats/GeometryVector/interface/Pi.h +++ b/DataFormats/GeometryVector/interface/Pi.h @@ -28,13 +28,13 @@ namespace Geom { - inline double pi() {return 3.141592653589793238;} - inline double twoPi() {return 2. *3.141592653589793238;} - inline double halfPi() {return 0.5*3.141592653589793238;} + inline constexpr double pi() {return 3.141592653589793238;} + inline constexpr double twoPi() {return 2. *3.141592653589793238;} + inline constexpr double halfPi() {return 0.5*3.141592653589793238;} - inline float fpi() {return 3.141592653589793238f;} - inline float ftwoPi() {return 2.f *3.141592653589793238f;} - inline float fhalfPi() {return 0.5f*3.141592653589793238f;} + inline constexpr float fpi() {return 3.141592653589793238f;} + inline constexpr float ftwoPi() {return 2.f *3.141592653589793238f;} + inline constexpr float fhalfPi() {return 0.5f*3.141592653589793238f;} } diff --git a/Geometry/CommonDetUnit/interface/GeomDet.h b/Geometry/CommonDetUnit/interface/GeomDet.h index f04384da60829..375ea2f51fa87 100644 --- a/Geometry/CommonDetUnit/interface/GeomDet.h +++ b/Geometry/CommonDetUnit/interface/GeomDet.h @@ -28,11 +28,13 @@ class SurfaceDeformation; class GeomDet { public: - typedef GeomDetEnumerators::SubDetector SubDetector; + using SubDetector = GeomDetEnumerators::SubDetector; + + + explicit GeomDet( Plane* plane): thePlane(plane) {} + explicit GeomDet( const ReferenceCountingPointer& plane) : thePlane(plane) {} - explicit GeomDet(Plane* plane); - explicit GeomDet(const ReferenceCountingPointer& plane); virtual ~GeomDet(); @@ -77,7 +79,7 @@ class GeomDet { DetId geographicalId() const { return m_detId; } /// Which subdetector - virtual SubDetector subDetector() const;; + virtual SubDetector subDetector() const; /// is a Unit virtual bool isLeaf() const { return components().empty();} @@ -93,10 +95,14 @@ class GeomDet { AlignmentPositionError const* alignmentPositionError() const { return theAlignmentPositionError;} - // specific unix index in a given subdetector (such as Tracker) + // specific unit index in a given subdetector (such as Tracker) int index() const { return m_index;} void setIndex(int i) { m_index=i;} + // specific geomDet index in a given subdetector (such as Tracker) + int gdetIndex() const { return m_gdetIndex;} + void setGdetIndex(int i) { m_gdetIndex=i;} + virtual const Topology& topology() const; @@ -105,7 +111,7 @@ class GeomDet { /// Return pointer to surface deformation. /// Defaults to "null" if not reimplemented in the derived classes. - virtual const SurfaceDeformation* surfaceDeformation() const { return 0; } + virtual const SurfaceDeformation* surfaceDeformation() const { return nullptr; } @@ -119,9 +125,10 @@ class GeomDet { ReferenceCountingPointer thePlane; DetId m_detId; - int m_index; + int m_index=-1; + int m_gdetIndex=-1; protected: - AlignmentPositionError* theAlignmentPositionError; + AlignmentPositionError* theAlignmentPositionError=nullptr; private: @@ -158,7 +165,7 @@ class GeomDet { }; -typedef GeomDet GeomDetUnit; +using GeomDetUnit = GeomDet; #endif diff --git a/Geometry/CommonDetUnit/src/GeomDet.cc b/Geometry/CommonDetUnit/src/GeomDet.cc index 3f4e6d2d6b671..d48a9545639fe 100644 --- a/Geometry/CommonDetUnit/src/GeomDet.cc +++ b/Geometry/CommonDetUnit/src/GeomDet.cc @@ -2,11 +2,6 @@ #include "Geometry/CommonDetUnit/interface/ModifiedSurfaceGenerator.h" #include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h" -GeomDet::GeomDet( Plane* plane): - thePlane(plane), m_index(-1), theAlignmentPositionError(nullptr) {} - -GeomDet::GeomDet( const ReferenceCountingPointer& plane) : - thePlane(plane), m_index(-1), theAlignmentPositionError(nullptr) {} GeomDet::~GeomDet() {delete theAlignmentPositionError;} diff --git a/Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h b/Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h index 170ceb7296001..075904fe0d723 100644 --- a/Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h +++ b/Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h @@ -45,14 +45,13 @@ class TrackerGeometry final : public TrackingGeometry { virtual ~TrackerGeometry() ; - - virtual const DetTypeContainer& detTypes() const; - virtual const DetUnitContainer& detUnits() const; - virtual const DetContainer& dets() const; - virtual const DetIdContainer& detUnitIds() const; - virtual const DetIdContainer& detIds() const; - virtual const TrackerGeomDet* idToDetUnit(DetId) const; - virtual const TrackerGeomDet* idToDet(DetId) const; + const DetTypeContainer& detTypes() const {return theDetTypes;} + const DetUnitContainer& detUnits() const {return theDetUnits;} + const DetContainer& dets() const {return theDets;} + const DetIdContainer& detUnitIds() const {return theDetUnitIds;} + const DetIdContainer& detIds() const { return theDetIds;} + const TrackerGeomDet* idToDetUnit(DetId) const; + const TrackerGeomDet* idToDet(DetId) const; const GeomDetEnumerators::SubDetector geomDetSubDetector(int subdet) const; unsigned int numberOfLayers(int subdet) const; diff --git a/Geometry/TrackerGeometryBuilder/src/TrackerGeometry.cc b/Geometry/TrackerGeometryBuilder/src/TrackerGeometry.cc index e0cbfd5552145..6a5d52b2647bb 100644 --- a/Geometry/TrackerGeometryBuilder/src/TrackerGeometry.cc +++ b/Geometry/TrackerGeometryBuilder/src/TrackerGeometry.cc @@ -156,6 +156,8 @@ void TrackerGeometry::addDetUnitId(DetId p){ } void TrackerGeometry::addDet(GeomDet const * p) { + // set index + const_cast(p)->setGdetIndex(theDets.size()); theDets.push_back(p); // add to vector theMap.insert(std::make_pair(p->geographicalId().rawId(),p)); DetId id(p->geographicalId()); @@ -189,17 +191,6 @@ void TrackerGeometry::addDetId(DetId p){ theDetIds.push_back(p); } -const TrackerGeometry::DetUnitContainer& -TrackerGeometry::detUnits() const -{ - return theDetUnits; -} - -const TrackerGeometry::DetContainer& -TrackerGeometry::dets() const -{ - return theDets; -} const TrackerGeometry::DetContainer& TrackerGeometry::detsPXB() const @@ -283,24 +274,6 @@ TrackerGeometry::isThere(GeomDetEnumerators::SubDetector subdet) const { return false; } -const TrackerGeometry::DetTypeContainer& -TrackerGeometry::detTypes() const -{ - return theDetTypes; -} - - -const TrackerGeometry::DetIdContainer& -TrackerGeometry::detUnitIds() const -{ - return theDetUnitIds; -} - -const TrackerGeometry::DetIdContainer& -TrackerGeometry::detIds() const -{ - return theDetIds; -} void TrackerGeometry::fillTestMap(const GeometricDet* gd) { std::string temp = gd->name(); diff --git a/RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h b/RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h index 7d4de2e2ffaf1..65383e8ceaed4 100644 --- a/RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h +++ b/RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h @@ -34,10 +34,10 @@ class ExceptionSafeStlPtrCol : public StlColType template RealType normalized_phi( RealType phi ) { - if (phi>CLHEP::pi) - { phi -= (2*CLHEP::pi) ; } - if (phi<-CLHEP::pi) - { phi += (2*CLHEP::pi) ; } + constexpr RealType pi(M_PI); + constexpr RealType pi2(2*M_PI); + if (phi>pi) { phi -= pi2 ; } + if (phi<-pi) { phi += pi2; } return phi ; } @@ -66,7 +66,7 @@ class EleRelPoint EleRelPoint( const GlobalPoint & p, const GlobalPoint & origin ) : relP_(p.x()-origin.x(),p.y()-origin.y(),p.z()-origin.z()) {} double eta() { return relP_.eta() ; } double phi() { return normalized_phi(relP_.phi()) ; } - double perp() { return sqrt(relP_.x()*relP_.x()+relP_.y()*relP_.y()) ; } + double perp() { return std::sqrt(relP_.x()*relP_.x()+relP_.y()*relP_.y()) ; } private : math::XYZVector relP_ ; } ; @@ -82,10 +82,10 @@ class EleRelPointPair EleRelPointPair( const math::XYZPoint & p1, const GlobalPoint & p2, const GlobalPoint & origin ) : relP1_(p1.x()-origin.x(),p1.y()-origin.y(),p1.z()-origin.z()), relP2_(p2.x()-origin.x(),p2.y()-origin.y(),p2.z()-origin.z()) {} EleRelPointPair( const GlobalPoint & p1, const math::XYZPoint & p2, const GlobalPoint & origin ) : relP1_(p1.x()-origin.x(),p1.y()-origin.y(),p1.z()-origin.z()), relP2_(p2.x()-origin.x(),p2.y()-origin.y(),p2.z()-origin.z()) {} EleRelPointPair( const GlobalPoint & p1, const GlobalPoint & p2, const GlobalPoint & origin ) : relP1_(p1.x()-origin.x(),p1.y()-origin.y(),p1.z()-origin.z()), relP2_(p2.x()-origin.x(),p2.y()-origin.y(),p2.z()-origin.z()) {} - double dEta() { return (relP1_.eta()-relP2_.eta()) ; } - double dPhi() { return normalized_phi(relP1_.phi()-relP2_.phi()) ; } - double dZ() { return (relP1_.z()-relP2_.z()) ; } - double dPerp() { return normalized_phi(relP1_.perp()-relP2_.perp()) ; } + auto dEta() { return (relP1_.eta()-relP2_.eta()) ; } + auto dPhi() { return normalized_phi(relP1_.barePhi()-relP2_.barePhi()) ; } + auto dZ() { return (relP1_.z()-relP2_.z()) ; } + auto dPerp() { return (relP1_.perp()-relP2_.perp()) ; } private : GlobalVector relP1_ ; GlobalVector relP2_ ; diff --git a/RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h b/RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h index 4d6e5251a0142..97de066daf4dc 100644 --- a/RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h +++ b/RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h @@ -231,7 +231,6 @@ class PixelHitMatcher bool searchInTIDTEC_ ; bool useRecoVertex_ ; - std::unordered_map mapTsos_fast_; std::unordered_map, TrajectoryStateOnSurface> mapTsos2_fast_; std::vector hit_gp_map_; } ; diff --git a/RecoEgamma/EgammaElectronAlgos/src/BarrelMeasurementEstimator.cc b/RecoEgamma/EgammaElectronAlgos/src/BarrelMeasurementEstimator.cc index 1254b1641334d..f51c234ecd296 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/BarrelMeasurementEstimator.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/BarrelMeasurementEstimator.cc @@ -18,7 +18,6 @@ #include "RecoEgamma/EgammaElectronAlgos/interface/BarrelMeasurementEstimator.h" #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h" -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h" #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h" #include "TrackingTools/DetLayers/interface/rangesIntersect.h" #include "TrackingTools/DetLayers/interface/PhiLess.h" @@ -81,20 +80,20 @@ std::pair BarrelMeasurementEstimator::estimate float zDiff = myZ -ts.z() ; float myZmax = theZMax; float myZmin = theZMin; - if(std::abs(myZ)<30. && myR>8.) + if( (std::abs(myZ)<30.f) & (myR>8.f) ) { - myZmax = 0.09; - myZmin = -0.09; + myZmax = 0.09f; + myZmin = -0.09f; } if( zDiff >= myZmax || zDiff <= myZmin ) return std::pair(false,0.); - float rhPhi = gp.phi() ; - float tsPhi = ts.phi(); + float rhPhi = gp.barePhi() ; + float tsPhi = ts.barePhi(); float phiDiff = normalized_phi(rhPhi-tsPhi) ; - if ( phiDiff < thePhiMax && phiDiff > thePhiMin ) + if ( (phiDiff < thePhiMax) & (phiDiff > thePhiMin) ) { return std::pair(true,1.) ; } else { return std::pair(false,0.) ; } @@ -108,13 +107,12 @@ bool BarrelMeasurementEstimator::estimate GlobalPoint trajPos(ts.globalParameters().position()); - GlobalDetRangeZPhi detRange(plane); Range trajZRange(trajPos.z() - std::abs(theZMin), trajPos.z() + std::abs(theZMax)); Range trajPhiRange(trajPos.phi() - std::abs(thePhiMin), trajPos.phi() + std::abs(thePhiMax)); - if(rangesIntersect(trajZRange, detRange.zRange()) && - rangesIntersect(trajPhiRange, detRange.phiRange(), PhiLess())) + if(rangesIntersect(trajZRange, plane.zSpan()) && + rangesIntersect(trajPhiRange, plane.phiSpan(), PhiLess())) { return true; } @@ -142,7 +140,7 @@ BarrelMeasurementEstimator::maximalLocalDisplacement if ( ts.hasError()) { LocalError le = ts.localError().positionError() ; - return Local2DVector( sqrt(le.xx())*nSigmaCut, sqrt(le.yy())*nSigmaCut) ; + return Local2DVector( std::sqrt(le.xx())*nSigmaCut, std::sqrt(le.yy())*nSigmaCut) ; } else return Local2DVector(99999,99999) ; } diff --git a/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc b/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc index bbc58c5c57ba8..19e8d37017b6a 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc @@ -47,6 +47,21 @@ #include #include + +#ifdef COUNT_ElectronSeeds +namespace { + struct Count { + long long s=0; + long long n=0; + ~Count() { std::cout << "ElectronSeeds res " << s<<'/'<("dynamicPhiRoad")), @@ -263,13 +278,19 @@ void ElectronSeedGenerator::run // Find the seeds recHits_.clear(); - LogDebug ("run") << "new cluster, calling seedsFromThisCluster"; + LogDebug ("ElectronSeedGenerator") << "new cluster, calling seedsFromThisCluster"; seedsFromThisCluster(sclRefs[i],hoe1s[i],hoe2s[i],out,tTopo); } - LogDebug ("run") << ": For event "< ForwardMeasurementEstimator::estimate( const TrajectoryStateOnSurface& ts, const TrackingRecHit& hit) const { @@ -38,14 +38,14 @@ std::pair ForwardMeasurementEstimator::estimate( const TrajectorySt float rMin = theRMin; float rMax = theRMax; float myZ = gp.z(); - if(std::abs(myZ)> 70. && std::abs(myZ)<170.) { + if( (std::abs(myZ)> 70.f) & (std::abs(myZ)<170.f)) { rMin = theRMinI; rMax = theRMaxI; } if( rDiff >= rMax || rDiff <= rMin ) return std::pair(false,0.); - float tsPhi = ts.globalParameters().position().phi(); - float rhPhi = gp.phi(); + float tsPhi = ts.globalParameters().position().barePhi(); + float rhPhi = gp.barePhi(); float myPhimin = thePhiMin; float myPhimax = thePhiMax; @@ -54,7 +54,7 @@ std::pair ForwardMeasurementEstimator::estimate( const TrajectorySt if (phiDiff > pi) phiDiff -= twopi; if (phiDiff < -pi) phiDiff += twopi; - if ( phiDiff < myPhimax && phiDiff > myPhimin ) { + if ( (phiDiff < myPhimax) & (phiDiff > myPhimin) ) { return std::pair(true,1.); } else { return std::pair(false,0.); @@ -75,15 +75,15 @@ std::pair ForwardMeasurementEstimator::estimate float rMin = theRMin; float rMax = theRMax; float myZ = gp.z(); - if(std::abs(myZ)> 70. && std::abs(myZ)<170.) { + if( (std::abs(myZ)> 70.f) & (std::abs(myZ)<170.f) ) { rMin = theRMinI; rMax = theRMaxI; } - if( rDiff >= rMax || rDiff <= rMin ) return std::pair(false,0.); + if( (rDiff >= rMax) | (rDiff <= rMin) ) return std::pair(false,0.); - float tsPhi = ts.phi(); - float rhPhi = gp.phi(); + float tsPhi = ts.barePhi(); + float rhPhi = gp.barePhi(); float myPhimin = thePhiMin; float myPhimax = thePhiMax; @@ -103,7 +103,6 @@ bool ForwardMeasurementEstimator::estimate typedef std::pair Range ; GlobalPoint trajPos(ts.globalParameters().position()); - GlobalDetRangeRPhi detRange(plane); float r1 = 0.; float r2 = 40.; @@ -111,8 +110,8 @@ bool ForwardMeasurementEstimator::estimate Range trajRRange(trajPos.perp() - r1, trajPos.perp() + r2); Range trajPhiRange(trajPos.phi() - std::abs(thePhiMin), trajPos.phi() + std::abs(thePhiMax)); - if(rangesIntersect(trajRRange, detRange.rRange()) && - rangesIntersect(trajPhiRange, detRange.phiRange(), PhiLess())) + if(rangesIntersect(trajRRange, plane.rSpan()) && + rangesIntersect(trajPhiRange, plane.phiSpan(), PhiLess())) { return true ; } else { return false ; } @@ -127,7 +126,7 @@ ForwardMeasurementEstimator::maximalLocalDisplacement if ( ts.hasError()) { LocalError le = ts.localError().positionError(); - return Local2DVector( sqrt(le.xx())*nSigmaCut, sqrt(le.yy())*nSigmaCut); + return Local2DVector( std::sqrt(le.xx())*nSigmaCut, std::sqrt(le.yy())*nSigmaCut); } else return Local2DVector(999999,999999) ; diff --git a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc index 1e7e80042ceb8..49a551810c4ce 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc @@ -436,7 +436,7 @@ void GsfElectronAlgo::ElectronData::computeCharge GlobalVector scvect(scpos-orig) ; GlobalPoint inntkpos = innTSOS.globalPosition() ; GlobalVector inntkvect = GlobalVector(inntkpos-orig) ; - float dPhiInnEle=normalized_phi(scvect.phi()-inntkvect.phi()) ; + float dPhiInnEle=normalized_phi(scvect.barePhi()-inntkvect.barePhi()) ; if(dPhiInnEle>0) info.scPixCharge = -1 ; else info.scPixCharge = 1 ; diff --git a/RecoEgamma/EgammaElectronAlgos/src/PixelHitMatcher.cc b/RecoEgamma/EgammaElectronAlgos/src/PixelHitMatcher.cc index 4251ecb39e933..dcb8e6063deee 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/PixelHitMatcher.cc +++ b/RecoEgamma/EgammaElectronAlgos/src/PixelHitMatcher.cc @@ -1,4 +1,3 @@ - #include "RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h" #include "RecoEgamma/EgammaElectronAlgos/interface/PixelMatchNextLayers.h" #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h" @@ -12,6 +11,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include +#include using namespace reco ; using namespace std ; @@ -94,8 +94,6 @@ vector PixelHitMatcher::predicted2Hits() float PixelHitMatcher::getVertex() { return vertex_ ; } -//CLHEP::Hep3Vector point_to_vector( const GlobalPoint & p ) -// { return CLHEP::Hep3Vector(p.x(),p.y(),p.z()) ; } std::vector PixelHitMatcher::compatibleSeeds @@ -106,52 +104,77 @@ PixelHitMatcher::compatibleSeeds typedef std::unordered_map, TrajectoryStateOnSurface> PosTsosAssoc; const int charge = int(fcharge) ; - const double xmeas_phi = xmeas.phi(); + // auto xmeas_phi = xmeas.barePhi(); + auto xmeas_r = xmeas.perp(); + + const float phicut = std::cos(2.5); + + FreeTrajectoryState fts = FTSFromVertexToPointFactory::get(*theMagField, xmeas, vprim, energy, charge); PerpendicularBoundPlaneBuilder bpb; TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum())); std::vector result ; - mapTsos_fast_.clear(); + //mapTsos_fast_.clear(); mapTsos2_fast_.clear(); - mapTsos_fast_.reserve(seeds->size()) ; + // mapTsos_fast_.reserve(seeds->size()) ; mapTsos2_fast_.reserve(seeds->size()) ; + // std::vector vTsos(theTrackerGeometry->dets().size()); + // TrajectoryStateOnSurface vTsos[theTrackerGeometry->dets().size()]; + + auto ndets = theTrackerGeometry->dets().size(); + + int iTsos[ndets]; + for ( auto & i : iTsos) i=-1; + std::vector vTsos; vTsos.reserve(seeds->size()); + for(const auto& seed : *seeds) { hit_gp_map_.clear(); if( seed.nHits() > 9 ) { edm::LogWarning("GsfElectronAlgo|UnexpectedSeed") <<"We cannot deal with seeds having more than 9 hits." ; continue; } + const TrajectorySeed::range& hits = seed.recHits(); // cache the global points + for( auto it = hits.first; it != hits.second; ++it ) { hit_gp_map_.emplace_back(it->globalPosition()); } - //iterate on the hits - for( auto it1 = hits.first; it1 != hits.second; ++it1 ) { + + //iterate on the hits + auto he = hits.second -1; + for( auto it1 = hits.first; it1 < he; ++it1 ) { if( !it1->isValid() ) continue; - const unsigned idx1 = std::distance(hits.first,it1); + auto idx1 = std::distance(hits.first,it1); const DetId id1 = it1->geographicalId(); - const GeomDet *geomdet1 = it1->det(); + const GeomDet *geomdet1 = it1->det(); + + auto ix1 = geomdet1->gdetIndex(); + + /* VI: this generates regression (other cut is just in phi). in my opinion it is safe and makes sense + auto away = geomdet1->position().basicVector().dot(xmeas.basicVector()) <0; + if (away) continue; + */ + const GlobalPoint& hit1Pos = hit_gp_map_[idx1]; + auto dt = hit1Pos.x()*xmeas.x()+hit1Pos.y()*xmeas.y(); + if (dt<0) continue; + if (dtsecond); - } else { - auto empl_result = - mapTsos_fast_.emplace(geomdet1,prop1stLayer->propagate(tsos,geomdet1->surface())); - tsos1 = &(empl_result.first->second); + if(iTsos[ix1]<0) { + iTsos[ix1] = vTsos.size(); + vTsos.push_back(prop1stLayer->propagate(tsos,geomdet1->surface())); } + auto tsos1 = &vTsos[iTsos[ix1]]; + if( !tsos1->isValid() ) continue; std::pair est = ( id1.subdetId() % 2 ? meas1stBLayer.estimate(vprim, *tsos1, hit1Pos) : meas1stFLayer.estimate(vprim, *tsos1, hit1Pos) ); if( !est.first ) continue; - if( std::abs(normalized_phi(hit1Pos.phi()-xmeas_phi))>2.5 ) continue; EleRelPointPair pp1(hit1Pos,tsos1->globalParameters().position(),vprim); const math::XYZPoint relHit1Pos(hit1Pos-vprim), relTSOSPos(tsos1->globalParameters().position() - vprim); const int subDet1 = id1.subdetId(); @@ -179,13 +202,13 @@ PixelHitMatcher::compatibleSeeds // now find the matching hit for( auto it2 = it1+1; it2 != hits.second; ++it2 ) { if( !it2->isValid() ) continue; - const unsigned idx2 = std::distance(hits.first,it2); + auto idx2 = std::distance(hits.first,it2); const DetId id2 = it2->geographicalId(); const GeomDet *geomdet2 = it2->det(); - const std::pair det_key(geomdet2,hit1Pos); + const std::pair det_key(geomdet2,hit1Pos); const TrajectoryStateOnSurface* tsos2; - PosTsosAssoc::iterator tsos2_itr = mapTsos2_fast_.find(det_key); - if( tsos2_itr != mapTsos2_fast_.end() ) { + auto tsos2_itr = mapTsos2_fast_.find(det_key); + if( tsos2_itr != mapTsos2_fast_.end() ) { tsos2 = &(tsos2_itr->second); } else { auto empl_result = @@ -209,9 +232,8 @@ PixelHitMatcher::compatibleSeeds }// outer loop on hits }// loop on seeds - mapTsos_fast_.clear() ; mapTsos2_fast_.clear() ; - + return result ; } @@ -256,7 +278,7 @@ PixelHitMatcher::compatibleHits LogDebug("") <<"[PixelHitMatcher::compatibleHits] nbr of hits compatible with extrapolation to first layer: " << pixelMeasurements.size(); for (aMeas m=pixelMeasurements.begin(); m!=pixelMeasurements.end(); m++){ if (m->recHit()->isValid()) { - float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ; + float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().barePhi()) ; if(std::abs(localDphi)>2.5)continue; CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(), m->forwardPredictedState().globalPosition().y(), @@ -286,7 +308,7 @@ PixelHitMatcher::compatibleHits for (aMeas m=pixel2Measurements.begin(); m!=pixel2Measurements.end(); m++){ if (m->recHit()->isValid()) { - float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ; + float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().barePhi()) ; if(std::abs(localDphi)>2.5)continue; CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(), m->forwardPredictedState().globalPosition().y(), @@ -327,7 +349,7 @@ PixelHitMatcher::compatibleHits for (aMeas m=pixelMeasurements.begin(); m!=pixelMeasurements.end(); m++){ if (m->recHit()->isValid()) { - float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()); + float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().barePhi()); if(std::abs(localDphi)>2.5)continue; CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(), m->forwardPredictedState().globalPosition().y(), @@ -348,7 +370,7 @@ PixelHitMatcher::compatibleHits for (aMeas m=pixel2Measurements.begin(); m!=pixel2Measurements.end(); m++){ if (m->recHit()->isValid()) { - float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ; + float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().barePhi()) ; if(std::abs(localDphi)>2.5)continue; CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(), m->forwardPredictedState().globalPosition().y(), @@ -413,7 +435,7 @@ PixelHitMatcher::compatibleHits if(!secondHit.measurementsInNextLayers().empty()){ for(unsigned int shit=0; shitglobalPosition().phi()) ; + float dphi = normalized_phi(pred1Meas[i].phi()-validMeasurements[i].recHit()->globalPosition().barePhi()) ; if (std::abs(dphi)<2.5) { ConstRecHitPointer pxrh = validMeasurements[i].recHit(); @@ -432,4 +454,3 @@ PixelHitMatcher::compatibleHits return result; } - diff --git a/RecoEgamma/EgammaPhotonAlgos/src/ConversionBarrelEstimator.cc b/RecoEgamma/EgammaPhotonAlgos/src/ConversionBarrelEstimator.cc index 054e503af756c..7a8fa396ba36c 100644 --- a/RecoEgamma/EgammaPhotonAlgos/src/ConversionBarrelEstimator.cc +++ b/RecoEgamma/EgammaPhotonAlgos/src/ConversionBarrelEstimator.cc @@ -5,7 +5,6 @@ #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" #include "TrackingTools/DetLayers/interface/PhiLess.h" #include "TrackingTools/DetLayers/interface/rangesIntersect.h" -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h" // zero value indicates incompatible ts - hit pair @@ -72,13 +71,12 @@ bool ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, // std::cout << " ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) " << std::endl; GlobalPoint trajPos(ts.globalParameters().position()); - GlobalDetRangeZPhi detRange(plane); Range trajZRange(trajPos.z() - 2.*theZRangeMax, trajPos.z() + 2.*theZRangeMax); Range trajPhiRange(trajPos.phi() + thePhiRangeMin, trajPos.phi() + thePhiRangeMax); - if(rangesIntersect(trajZRange, detRange.zRange()) && - rangesIntersect(trajPhiRange, detRange.phiRange(), PhiLess())) { + if(rangesIntersect(trajZRange, plane.zSpan()) && + rangesIntersect(trajPhiRange, plane.phiSpan(), PhiLess())) { // std::cout << " ConversionBarrelEstimator::estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane) IN RANGE " << std::endl; return true; diff --git a/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeRPhi.h b/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeRPhi.h deleted file mode 100644 index 25475a0e196e7..0000000000000 --- a/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeRPhi.h +++ /dev/null @@ -1,25 +0,0 @@ -#ifndef RecoTracker_TkTrackingRegions_GlobalDetRangeRPhi_H -#define RecoTracker_TkTrackingRegions_GlobalDetRangeRPhi_H - -#include - -class Plane; - -/** Keep R and Phi range for detunit */ - -class GlobalDetRangeRPhi { -public: - - typedef std::pair Range; - - GlobalDetRangeRPhi( const Plane& det); - - Range rRange() const { return theRRange;} - Range phiRange() const { return thePhiRange;} - -private: - Range theRRange; - Range thePhiRange; -}; - -#endif diff --git a/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h b/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h deleted file mode 100644 index 9dca8baeec95e..0000000000000 --- a/RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h +++ /dev/null @@ -1,26 +0,0 @@ -#ifndef RecoTracker_TkTrackingRegions_GlobalDetRangeZPhi_H -#define RecoTracker_TkTrackingRegions_GlobalDetRangeZPhi_H - -#include - -class Plane; - -/** Implementation class for PhiZMeasurementEstimator etc. - */ - -class GlobalDetRangeZPhi { -public: - - typedef std::pair Range; - - GlobalDetRangeZPhi( const Plane& det); - - Range zRange() const { return theZRange;} - Range phiRange() const { return thePhiRange;} - -private: - Range theZRange; - Range thePhiRange; -}; - -#endif diff --git a/RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h b/RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h index e321a77414c71..ee8f04aeefe0a 100644 --- a/RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h +++ b/RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h @@ -5,8 +5,6 @@ #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h" #include "RecoTracker/TkTrackingRegions/interface/TkTrackingRegionsMargin.h" -#include "FWCore/Utilities/interface/GCC11Compatibility.h" - class HitRZCompatibility { public: // only three algos are implemented.. diff --git a/RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h b/RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h index 6c10cea9cc8db..5769755e46ba3 100644 --- a/RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h +++ b/RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h @@ -21,7 +21,6 @@ class OuterHitPhiPrediction; -class OuterEstimator; class BarrelDetLayer; class ForwardDetLayer; class MeasurementTrackerEvent; @@ -139,9 +138,10 @@ class RectangularEtaPhiTrackingRegion final : public TrackingRegion { UseMeasurementTracker whereToUseMeasurementTracker = UseMeasurementTracker::kNever, bool precise = true, const MeasurementTrackerEvent *measurementTracker = nullptr, - bool etaPhiRegion=false) + bool etaPhiRegion=false, bool useMS=true) : TrackingRegionBase( dir, vertexPos, invPtRange, rVertex, zVertex), - thePhiMargin( phiMargin), theMeasurementTrackerUsage(whereToUseMeasurementTracker), thePrecise(precise),theUseEtaPhi(etaPhiRegion), + thePhiMargin( phiMargin), theMeasurementTrackerUsage(whereToUseMeasurementTracker), + thePrecise(precise),theUseMS(useMS),theUseEtaPhi(etaPhiRegion), theMeasurementTracker(measurementTracker) { initEtaRange(dir, etaMargin); } @@ -163,10 +163,10 @@ class RectangularEtaPhiTrackingRegion final : public TrackingRegion { virtual HitRZCompatibility * checkRZ(const DetLayer* layer, const Hit & outerHit, - const edm::EventSetup& iSetup, - const DetLayer* outerlayer=0, + const edm::EventSetup&iSetup, + const DetLayer* outerlayer=nullptr, float lr=0, float gz=0, float dr=0, float dz=0) const - { return checkRZOld(layer,outerHit->hit(),iSetup); } + { return checkRZOld(layer,outerHit,iSetup, outerlayer); } virtual RectangularEtaPhiTrackingRegion* clone() const { return new RectangularEtaPhiTrackingRegion(*this); @@ -178,11 +178,12 @@ class RectangularEtaPhiTrackingRegion final : public TrackingRegion { private: HitRZCompatibility* checkRZOld( const DetLayer* layer, - const TrackingRecHit* outerHit, - const edm::EventSetup& iSetup) const; + const Hit & outerHit, + const edm::EventSetup&iSetup, + const DetLayer* outerlayer) const; - std::unique_ptr estimator(const BarrelDetLayer* layer,const edm::EventSetup& iSetup) const dso_internal; - std::unique_ptr estimator(const ForwardDetLayer* layer,const edm::EventSetup& iSetup) const dso_internal; + std::unique_ptr estimator(const BarrelDetLayer* layer,const edm::EventSetup& iSetup) const dso_internal; + std::unique_ptr estimator(const ForwardDetLayer* layer,const edm::EventSetup& iSetup) const dso_internal; OuterHitPhiPrediction phiWindow(const edm::EventSetup& iSetup) const dso_internal; HitRZConstraint rzConstraint() const dso_internal; @@ -195,10 +196,11 @@ class RectangularEtaPhiTrackingRegion final : public TrackingRegion { Range theLambdaRange; Margin thePhiMargin; float theMeanLambda; - const UseMeasurementTracker theMeasurementTrackerUsage; - bool thePrecise; - bool theUseEtaPhi; - const MeasurementTrackerEvent *theMeasurementTracker; + const UseMeasurementTracker theMeasurementTrackerUsage = UseMeasurementTracker::kNever; + bool thePrecise=false; + bool theUseMS=false; + bool theUseEtaPhi=false; + const MeasurementTrackerEvent *theMeasurementTracker = nullptr; diff --git a/RecoTracker/TkTrackingRegions/src/GlobalDetRangeRPhi.cc b/RecoTracker/TkTrackingRegions/src/GlobalDetRangeRPhi.cc deleted file mode 100644 index e817386ae4ed9..0000000000000 --- a/RecoTracker/TkTrackingRegions/src/GlobalDetRangeRPhi.cc +++ /dev/null @@ -1,42 +0,0 @@ -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeRPhi.h" -#include "DataFormats/GeometrySurface/interface/BoundPlane.h" -#include "DataFormats/GeometrySurface/interface/Bounds.h" -#include "TrackingTools/DetLayers/interface/PhiLess.h" -#include -#include -#include - -using namespace std; - -GlobalDetRangeRPhi::GlobalDetRangeRPhi( const BoundPlane& plane) { - - float dx = plane.bounds().width()/2.; - float dy = plane.bounds().length()/2.; - - std::vector corners(4); - corners[0] = plane.toGlobal( LocalPoint( -dx, -dy, 0)); - corners[1] = plane.toGlobal( LocalPoint( -dx, dy, 0)); - corners[2] = plane.toGlobal( LocalPoint( dx, -dy, 0)); - corners[3] = plane.toGlobal( LocalPoint( dx, dy, 0)); - - float phimin = corners[0].phi(); float phimax = phimin; - float rmin = corners[0].perp(); float rmax = rmin; - for ( int i=1; i<4; i++) { - float phi = corners[i].phi(); - if ( PhiLess()( phi, phimin)) phimin = phi; - if ( PhiLess()( phimax, phi)) phimax = phi; - -// if ( phi < phimin) phimin = phi; -// if ( phi > phimax) phimax = phi; - - float r = corners[i].perp(); - if ( r < rmin) rmin = r; - if ( r > rmax) rmax = r; - } - - theRRange.first = rmin; - theRRange.second = rmax; - thePhiRange.first = phimin; - thePhiRange.second = phimax; - -} diff --git a/RecoTracker/TkTrackingRegions/src/GlobalDetRangeZPhi.cc b/RecoTracker/TkTrackingRegions/src/GlobalDetRangeZPhi.cc deleted file mode 100644 index 510a5f61c0c9a..0000000000000 --- a/RecoTracker/TkTrackingRegions/src/GlobalDetRangeZPhi.cc +++ /dev/null @@ -1,39 +0,0 @@ -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h" -#include "DataFormats/GeometrySurface/interface/BoundPlane.h" -#include "DataFormats/GeometrySurface/interface/Bounds.h" -#include "TrackingTools/DetLayers/interface/PhiLess.h" -#include -#include -#include - -using namespace std; - -GlobalDetRangeZPhi::GlobalDetRangeZPhi( const BoundPlane& plane) { - - float dx = plane.bounds().width()/2.; - float dy = plane.bounds().length()/2.; - - vector corners(4); - corners[0] = plane.toGlobal( LocalPoint( -dx, -dy, 0)); - corners[1] = plane.toGlobal( LocalPoint( -dx, dy, 0)); - corners[2] = plane.toGlobal( LocalPoint( dx, -dy, 0)); - corners[3] = plane.toGlobal( LocalPoint( dx, dy, 0)); - - float phimin = corners[0].phi(); float phimax = phimin; - float zmin = corners[0].z(); float zmax = zmin; - for ( int i=1; i<4; i++) { - float phi = corners[i].phi(); - if ( PhiLess()( phi, phimin)) phimin = phi; - if ( PhiLess()( phimax, phi)) phimax = phi; - - float z = corners[i].z(); - if ( z < zmin) zmin = z; - if ( z > zmax) zmax = z; - } - - theZRange.first = zmin; - theZRange.second = zmax; - thePhiRange.first = phimin; - thePhiRange.second = phimax; - -} diff --git a/RecoTracker/TkTrackingRegions/src/OuterDetCompatibility.cc b/RecoTracker/TkTrackingRegions/src/OuterDetCompatibility.cc index 7834ec9aa5190..471ba7ccf2820 100644 --- a/RecoTracker/TkTrackingRegions/src/OuterDetCompatibility.cc +++ b/RecoTracker/TkTrackingRegions/src/OuterDetCompatibility.cc @@ -1,6 +1,4 @@ #include "OuterDetCompatibility.h" -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeZPhi.h" -#include "RecoTracker/TkTrackingRegions/interface/GlobalDetRangeRPhi.h" #include "TrackingTools/DetLayers/interface/PhiLess.h" #include "TrackingTools/DetLayers/interface/rangesIntersect.h" @@ -9,13 +7,11 @@ using namespace std; bool OuterDetCompatibility::operator() (const BoundPlane& plane) const { if (barrel) { - GlobalDetRangeZPhi detRange(plane); - if (!checkPhi(detRange.phiRange())) return 0; - if (!checkZ(detRange.zRange())) return 0; + if (!checkPhi(plane.phiSpan())) return false; + if (!checkZ(plane.zSpan())) return false; } else { - GlobalDetRangeRPhi detRange(plane); - if (!checkPhi(detRange.phiRange())) return 0; - if (!checkR(detRange.rRange())) return 0; + if (!checkPhi(plane.phiSpan())) return false; + if (!checkR(plane.rSpan())) return false; } return 1; } diff --git a/RecoTracker/TkTrackingRegions/src/OuterEstimator.h b/RecoTracker/TkTrackingRegions/src/OuterEstimator.h index 2398c62cf6575..ab3766120b80f 100644 --- a/RecoTracker/TkTrackingRegions/src/OuterEstimator.h +++ b/RecoTracker/TkTrackingRegions/src/OuterEstimator.h @@ -21,50 +21,56 @@ #include "FWCore/Utilities/interface/Visibility.h" +template class dso_internal OuterEstimator final : public MeasurementEstimator { public: + + using OuterHitCompat = OuterHitCompatibility; + OuterEstimator( const OuterDetCompatibility & detCompatibility, - const OuterHitCompatibility & hitCompatibility, + const OuterHitCompat & hitCompatibility, const edm::EventSetup& iSetup) : theDetCompatibility(detCompatibility), theHitCompatibility (hitCompatibility) { } - virtual ~OuterEstimator(){} - virtual std::pair estimate( + + ~OuterEstimator(){} + + std::pair estimate( const TrajectoryStateOnSurface& ts, const TrackingRecHit& hit) - const { + const override { return theHitCompatibility(hit) ? std::make_pair(true,1.) : std::make_pair(false,0.) ; } - virtual bool estimate( + bool estimate( const TrajectoryStateOnSurface& ts, const BoundPlane& plane -) const { + ) const override { return theDetCompatibility(plane); } GlobalPoint center() { return theDetCompatibility.center(); } - virtual OuterEstimator* clone() const { + OuterEstimator* clone() const override { return new OuterEstimator(*this); } - virtual MeasurementEstimator::Local2DVector maximalLocalDisplacement( - const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const { + MeasurementEstimator::Local2DVector maximalLocalDisplacement( + const TrajectoryStateOnSurface& ts, const BoundPlane& plane) const override { return theDetCompatibility.maximalLocalDisplacement( ts.globalPosition(),plane); } const OuterDetCompatibility & detCompatibility() const {return theDetCompatibility; } - const OuterHitCompatibility & hitCompatibility() const + const OuterHitCompat & hitCompatibility() const {return theHitCompatibility; } private: OuterDetCompatibility theDetCompatibility; - OuterHitCompatibility theHitCompatibility; + OuterHitCompat theHitCompatibility; }; #endif diff --git a/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.cc b/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.cc deleted file mode 100644 index 8bc327a37ce8f..0000000000000 --- a/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.cc +++ /dev/null @@ -1,25 +0,0 @@ -#include "OuterHitCompatibility.h" -#include "TrackingTools/DetLayers/interface/PhiLess.h" -#include "DataFormats/Common/interface/Handle.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" - -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDetType.h" -#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" - -bool OuterHitCompatibility::operator() ( const TrackingRecHit & hit) const -{ - GlobalPoint hitPos = hit.globalPosition(); - float hitR = hitPos.perp(); - float hitPhi = hitPos.phi(); - - if ( !checkPhi(hitPhi, hitR) ) return false; - - float hitZ = hitPos.z(); - if ( !(*theRZCompatibility)(hitR,hitZ) ) return false; - - return true; -} - diff --git a/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.h b/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.h index efd8bad715bcc..ca9aa540ebd62 100644 --- a/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.h +++ b/RecoTracker/TkTrackingRegions/src/OuterHitCompatibility.h @@ -15,40 +15,43 @@ #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" #include "TrackingTools/DetLayers/interface/PhiLess.h" +#include "DataFormats/Math/interface/approx_atan2.h" + #include "FWCore/Utilities/interface/Visibility.h" +template class dso_internal OuterHitCompatibility { public: OuterHitCompatibility( const OuterHitPhiPrediction & phiPrediction, - const HitRZCompatibility & rzCompatibility) - : thePhiPrediction(phiPrediction) - { theRZCompatibility = rzCompatibility.clone(); } + const Algo & rzCompatibility) + : thePhiPrediction(phiPrediction), + theRZCompatibility(rzCompatibility) {} + + bool operator() (const TrackingRecHit & hit) const { + + auto hitPos = hit.globalPosition(); + auto hitR = hitPos.perp(); - OuterHitCompatibility(const OuterHitCompatibility & ohc) - : thePhiPrediction(ohc.thePhiPrediction) - { theRZCompatibility = ohc.theRZCompatibility->clone(); } + auto hitZ = hitPos.z(); + if ( !theRZCompatibility(hitR,hitZ) ) return false; - ~OuterHitCompatibility() - { delete theRZCompatibility; } + auto hitPhi = unsafe_atan2f<9>(hitPos.y(),hitPos.x()); + return checkPhi(hitPhi, hitR); + } - bool operator() (const TrackingRecHit & hit) const; bool checkPhi(float phi, float r) const { - OuterHitPhiPrediction::Range hitPhiRange = thePhiPrediction(r); + auto hitPhiRange = thePhiPrediction(r); PhiLess less; bool phiOK = less(hitPhiRange.min(),phi) && less(phi,hitPhiRange.max()); return phiOK; } - OuterHitCompatibility* clone() const { - return new OuterHitCompatibility(*this); - } - -protected: - const HitRZCompatibility * theRZCompatibility; +private: OuterHitPhiPrediction thePhiPrediction; + Algo theRZCompatibility; }; #endif diff --git a/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.cc b/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.cc index 3431be9df6e6c..fe55299204d4a 100644 --- a/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.cc +++ b/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.cc @@ -1,24 +1,37 @@ -#include #include "OuterHitPhiPrediction.h" +#include "DataFormats/Math/interface/approx_asin.h" OuterHitPhiPrediction::Range - OuterHitPhiPrediction::operator()(float radius) const + OuterHitPhiPrediction::sym(float radius) const +{ + auto arc = radius*theCurvature.max()*0.5f + theOriginRBound/radius; + + auto Phi_r = unsafe_asin07<5>(arc); + return Range( thePhiAtVertex.min() - Phi_r - theTolerance, + thePhiAtVertex.max() + Phi_r + theTolerance); +} + + +// in case somebody comes with a RELEVANT use case... +OuterHitPhiPrediction::Range + OuterHitPhiPrediction::asym(float radius) const { auto invr = 1.f/radius; - if( std::max(std::abs(theCurvature.min()), std::abs(theCurvature.max())) > invr) - return Range(-M_PI,M_PI); - + if( std::max(std::abs(theCurvature.min()), std::abs(theCurvature.max())) > invr) + return Range(-M_PI,M_PI); + float Phi_r = std::asin(radius*theCurvature.max()*0.5f + theOriginRBound*invr); if (theCurvature.max() == -theCurvature.min()) - return Range( thePhiAtVertex.min() - Phi_r - theTolerance.left(), - thePhiAtVertex.max() + Phi_r + theTolerance.right()); + return Range( thePhiAtVertex.min() - Phi_r - theTolerance, + thePhiAtVertex.max() + Phi_r + theTolerance); float curv0 = theCurvature.mean(); float Phi_0 = std::asin(radius*curv0*0.5f); float Phi_m = std::asin(radius*theCurvature.min()*0.5f-theOriginRBound*invr); - return Range( thePhiAtVertex.min() + Phi_0 + Phi_m - theTolerance.left(), - thePhiAtVertex.max() + Phi_0 + Phi_r + theTolerance.right()); + return Range( thePhiAtVertex.min() + Phi_0 + Phi_m - theTolerance, + thePhiAtVertex.max() + Phi_0 + Phi_r + theTolerance); } + diff --git a/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.h b/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.h index c7689e89ff297..5129422aa1346 100644 --- a/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.h +++ b/RecoTracker/TkTrackingRegions/src/OuterHitPhiPrediction.h @@ -4,32 +4,40 @@ /** predicts phi range at a given radius r */ #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h" -#include "RecoTracker/TkTrackingRegions/interface/TkTrackingRegionsMargin.h" #include "FWCore/Utilities/interface/Visibility.h" +#include + class dso_internal OuterHitPhiPrediction { public: - typedef PixelRecoRange Range; - typedef TkTrackingRegionsMargin Margin; + using Range= PixelRecoRange; OuterHitPhiPrediction( const Range & phiAtVertex, const Range & curvature, - float originRBound, - const Margin & tolerance = Margin(0,0)) + float originRBound) : thePhiAtVertex(phiAtVertex), theCurvature(curvature), - theOriginRBound (originRBound), theTolerance(tolerance) { } + theOriginRBound (originRBound) { + assert(theCurvature.max()>0); + assert(theCurvature.max() == -theCurvature.min()); + } + + void setTolerance(float tolerance) { theTolerance = tolerance; } - void setTolerance(const Margin & tolerance) { theTolerance = tolerance; } - Range operator()(float radius) const; + Range operator()(float radius) const { return sym(radius);} private: + + Range sym(float radius) const; + Range asym(float radius) const; + + Range thePhiAtVertex; Range theCurvature; float theOriginRBound; - Margin theTolerance; + float theTolerance = 0.f; }; #endif diff --git a/RecoTracker/TkTrackingRegions/src/RectangularEtaPhiTrackingRegion.cc b/RecoTracker/TkTrackingRegions/src/RectangularEtaPhiTrackingRegion.cc index 0af0077a0f0c6..0e2c3ba251fd9 100644 --- a/RecoTracker/TkTrackingRegions/src/RectangularEtaPhiTrackingRegion.cc +++ b/RecoTracker/TkTrackingRegions/src/RectangularEtaPhiTrackingRegion.cc @@ -67,15 +67,12 @@ void RectangularEtaPhiTrackingRegion:: initEtaRange( const GlobalVector & dir, c } HitRZCompatibility* RectangularEtaPhiTrackingRegion:: -checkRZOld(const DetLayer* layer, const TrackingRecHit *outerHit,const edm::EventSetup& iSetup) const +checkRZOld(const DetLayer* layer, const Hit & outerHit, const edm::EventSetup&iSetup, const DetLayer* outerlayer) const { - edm::ESHandle tracker; - iSetup.get().get(tracker); bool isBarrel = (layer->location() == GeomDetEnumerators::barrel); - GlobalPoint ohit = tracker->idToDet(outerHit->geographicalId())->surface().toGlobal(outerHit->localPosition()); - float outerred_r = sqrt( sqr(ohit.x()-origin().x())+sqr(ohit.y()-origin().y()) ); - //PixelRecoPointRZ outer(ohit.perp(), ohit.z()); + GlobalPoint ohit = outerHit->globalPosition(); + float outerred_r = std::sqrt( sqr(ohit.x()-origin().x())+sqr(ohit.y()-origin().y()) ); PixelRecoPointRZ outer(outerred_r, ohit.z()); float zMinOrigin = origin().z() - originZBound(); @@ -92,61 +89,63 @@ checkRZOld(const DetLayer* layer, const TrackingRecHit *outerHit,const edm::Even float cotLeft = std::min(vcotMax, theLambdaRange.max()); return new HitEtaCheck( isBarrel, outer, cotLeft, cotRight); } - float hitZErr = 0.; - float hitRErr = 0.; - - PixelRecoPointRZ outerL, outerR; - if (layer->location() == GeomDetEnumerators::barrel) { - outerL = PixelRecoPointRZ(outer.r(), outer.z()-hitZErr); - outerR = PixelRecoPointRZ(outer.r(), outer.z()+hitZErr); - } else if (outer.z() > 0) { - outerL = PixelRecoPointRZ(outer.r()+hitRErr, outer.z()); - outerR = PixelRecoPointRZ(outer.r()-hitRErr, outer.z()); - } else { - outerL = PixelRecoPointRZ(outer.r()-hitRErr, outer.z()); - outerR = PixelRecoPointRZ(outer.r()+hitRErr, outer.z()); - } + + + + float outerZscatt=0; + float innerScatt =0; //CHECK - MultipleScatteringParametrisation oSigma(layer,iSetup); - float cotThetaOuter = theMeanLambda; - float sinThetaOuter = 1/std::sqrt(1+sqr(cotThetaOuter)); - float outerZscatt = 3.f*oSigma(ptMin(),cotThetaOuter) / sinThetaOuter; + if (theUseMS) { + MultipleScatteringParametrisation oSigma(layer,iSetup); + float cotThetaOuter = theMeanLambda; + float sinThetaOuterInv = std::sqrt(1.f+sqr(cotThetaOuter)); + outerZscatt = 3.f*oSigma(ptMin(),cotThetaOuter)*sinThetaOuterInv; + } - PixelRecoLineRZ boundL(outerL, theLambdaRange.max()); - PixelRecoLineRZ boundR(outerR, theLambdaRange.min()); + PixelRecoLineRZ boundL(outer, theLambdaRange.max()); + PixelRecoLineRZ boundR(outer, theLambdaRange.min()); float zMinLine = boundL.zAtR(0.)-outerZscatt; float zMaxLine = boundR.zAtR(0.)+outerZscatt; PixelRecoPointRZ vtxL(0.,max(zMinLine, zMinOrigin)); PixelRecoPointRZ vtxR(0.,min(zMaxLine, zMaxOrigin)); PixelRecoPointRZ vtxMean(0.,(vtxL.z()+vtxR.z())*0.5f); //CHECK - MultipleScatteringParametrisation iSigma(layer,iSetup); - float innerScatt = 3.f * iSigma(ptMin(),vtxMean, outer); - - SimpleLineRZ leftLine( vtxL, outerL); - SimpleLineRZ rightLine( vtxR, outerR); + + if (theUseMS) { + MultipleScatteringParametrisation iSigma(layer,iSetup); + + innerScatt = 3.f * ( outerlayer ? + iSigma( ptMin(), vtxMean, outer, outerlayer->seqNum()) + : iSigma( ptMin(), vtxMean, outer) ) ; + + // innerScatt = 3.f *iSigma( ptMin(), vtxMean, outer); + } + + SimpleLineRZ leftLine( vtxL, outer); + SimpleLineRZ rightLine(vtxR, outer); HitRZConstraint rzConstraint(leftLine, rightLine); - float cotTheta = std::abs(leftLine.cotLine()+rightLine.cotLine())*0.5f; + auto cotTheta = std::abs(leftLine.cotLine()+rightLine.cotLine())*0.5f; -// float bendR = longitudinalBendingCorrection(outer.r(),ptMin()); // std::cout << "RectangularEtaPhiTrackingRegion " << outer.r()<<','<< outer.z() << " " << innerScatt << " " << cotTheta << " " << hitZErr << std::endl; if (isBarrel) { - float sinTheta = 1/std::sqrt(1+sqr(cotTheta)); - float corrZ = innerScatt/sinTheta + hitZErr; - return new HitZCheck(rzConstraint, HitZCheck::Margin(corrZ,corrZ)); + auto sinThetaInv = std::sqrt(1.f+sqr(cotTheta)); + auto corr = innerScatt*sinThetaInv; + return new HitZCheck(rzConstraint, HitZCheck::Margin(corr,corr)); } else { - float cosTheta = 1/std::sqrt(1+sqr(1/cotTheta)); - float corrR = innerScatt/cosTheta + hitRErr; - return new HitRCheck( rzConstraint, HitRCheck::Margin(corrR,corrR)); + auto cosThetaInv = std::sqrt(1.f+sqr(1.f/cotTheta)); + auto corr = innerScatt*cosThetaInv; + return new HitRCheck( rzConstraint, HitRCheck::Margin(corr,corr)); } } -std::unique_ptr +std::unique_ptr RectangularEtaPhiTrackingRegion::estimator(const BarrelDetLayer* layer,const edm::EventSetup& iSetup) const { + + using Algo = HitZCheck; // det dimensions float halfLength = 0.5f*layer->surface().bounds().length(); @@ -162,7 +161,7 @@ std::unique_ptr HitZCheck zPrediction(rzConstraint()); Range hitZWindow = zPrediction.range(detRWindow.min()). intersection(detZWindow); - if (hitZWindow.empty()) return 0; + if (hitZWindow.empty()) return nullptr; // phi prediction OuterHitPhiPrediction phiPrediction = phiWindow(iSetup); @@ -172,18 +171,19 @@ std::unique_ptr // OuterHitPhiPrediction::Range phiRange; if (thePrecise) { - float cotTheta = (hitZWindow.mean()-origin().z()) / radius; - float sinTheta = 1/std::sqrt(1+sqr(cotTheta)); + auto invR = 1.f/ radius; + auto cotTheta = (hitZWindow.mean()-origin().z()) * invR; + auto sinThetaInv = std::sqrt(1.f+sqr(cotTheta)); MultipleScatteringParametrisation msSigma(layer,iSetup); - float scatt = 3.f * msSigma(ptMin(), cotTheta); - float bendR = longitudinalBendingCorrection(radius,ptMin(),iSetup); + auto scatt = 3.f * msSigma(ptMin(), cotTheta); + auto bendR = longitudinalBendingCorrection(radius,ptMin(),iSetup); float hitErrRPhi = 0.; float hitErrZ = 0.; - float corrPhi = (scatt+ hitErrRPhi)/radius; - float corrZ = scatt/sinTheta + bendR*std::abs(cotTheta) + hitErrZ; + float corrPhi = (scatt+ hitErrRPhi)*invR; + float corrZ = scatt*sinThetaInv + bendR*std::abs(cotTheta) + hitErrZ; - phiPrediction.setTolerance(OuterHitPhiPrediction::Margin(corrPhi,corrPhi)); + phiPrediction.setTolerance(corrPhi); zPrediction.setTolerance(HitZCheck::Margin(corrZ,corrZ)); // @@ -200,16 +200,16 @@ std::unique_ptr phiRange = phiPrediction(detRWindow.mean()); } - return std::make_unique( + return std::make_unique>( OuterDetCompatibility( layer, phiRange, detRWindow, hitZWindow), - OuterHitCompatibility( phiPrediction, zPrediction ), + OuterHitCompatibility( phiPrediction, zPrediction ), iSetup); } -std::unique_ptr +std::unique_ptr RectangularEtaPhiTrackingRegion::estimator(const ForwardDetLayer* layer,const edm::EventSetup& iSetup) const { - + using Algo=HitRCheck; // det dimensions, ranges float halfThickness = 0.5f*layer->surface().bounds().thickness(); float zLayer = layer->position().z() ; @@ -231,16 +231,16 @@ RectangularEtaPhiTrackingRegion::estimator(const ForwardDetLayer* layer,const ed // if (thePrecise) { float cotTheta = (detZWindow.mean()-origin().z())/hitRWindow.mean(); - float cosTheta = cotTheta/std::sqrt(1+sqr(cotTheta)); + float cosThetaInv = std::sqrt(1+sqr(cotTheta))/cotTheta; MultipleScatteringParametrisation msSigma(layer,iSetup); float scatt = 3.f * msSigma(ptMin(),cotTheta); float bendR = longitudinalBendingCorrection(hitRWindow.max(),ptMin(),iSetup); float hitErrRPhi = 0.; float hitErrR = 0.; float corrPhi = (scatt+hitErrRPhi)/detRWindow.min(); - float corrR = scatt/std::abs(cosTheta) + bendR + hitErrR; + float corrR = scatt*std::abs(cosThetaInv) + bendR + hitErrR; - phiPrediction.setTolerance(OuterHitPhiPrediction::Margin(corrPhi,corrPhi)); + phiPrediction.setTolerance(corrPhi); rPrediction.setTolerance(HitRCheck::Margin(corrR,corrR)); // @@ -257,9 +257,9 @@ RectangularEtaPhiTrackingRegion::estimator(const ForwardDetLayer* layer,const ed hitRWindow = Range(w1.min(),w2.max()).intersection(detRWindow); } - return std::make_unique( + return std::make_unique>( OuterDetCompatibility( layer, phiRange, hitRWindow, detZWindow), - OuterHitCompatibility( phiPrediction, rPrediction),iSetup ); + OuterHitCompatibility( phiPrediction, rPrediction),iSetup ); } @@ -267,7 +267,7 @@ RectangularEtaPhiTrackingRegion::estimator(const ForwardDetLayer* layer,const ed OuterHitPhiPrediction RectangularEtaPhiTrackingRegion::phiWindow(const edm::EventSetup& iSetup) const { - float phi0 = phiDirection(); + auto phi0 = phiDirection(); return OuterHitPhiPrediction( OuterHitPhiPrediction::Range( phi0-thePhiMargin.left(), phi0+thePhiMargin.right()), @@ -308,7 +308,6 @@ TrackingRegion::Hits RectangularEtaPhiTrackingRegion::hits( //ESTIMATOR const DetLayer * detLayer = layer.detLayer(); - std::unique_ptr est; bool measurementMethod = false; if(theMeasurementTrackerUsage == UseMeasurementTracker::kAlways) measurementMethod = true; @@ -323,6 +322,7 @@ TrackingRegion::Hits RectangularEtaPhiTrackingRegion::hits( const GlobalPoint vtx = origin(); GlobalVector dir = direction(); + std::unique_ptr est; if ((GeomDetEnumerators::isTrackerPixel(detLayer->subDetector()) && GeomDetEnumerators::isBarrel(detLayer->subDetector())) || (!theUseEtaPhi && detLayer->location() == GeomDetEnumerators::barrel)) { const BarrelDetLayer& bl = dynamic_cast(*detLayer); @@ -386,24 +386,36 @@ TrackingRegion::Hits RectangularEtaPhiTrackingRegion::hits( } else { // - // temporary solution + // temporary solution (actually heavily used for Pixels....) // if (detLayer->location() == GeomDetEnumerators::barrel) { const BarrelDetLayer& bl = dynamic_cast(*detLayer); - est = estimator(&bl,es); + auto est = estimator(&bl,es); + if (!est) return result; + using Algo = HitZCheck; + auto const & hitComp = (reinterpret_cast const&>(*est)).hitCompatibility(); + auto layerHits = layer.hits(); + result.reserve(layerHits.size()); + for (auto && ih : layerHits) { + if ( hitComp(*ih) ) + result.emplace_back( std::move(ih) ); + } + } else { const ForwardDetLayer& fl = dynamic_cast(*detLayer); - est = estimator(&fl,es); - } - if (!est) return result; - - auto layerHits = layer.hits(); - result.reserve(layerHits.size()); - for (auto && ih : layerHits) { - if ( est->hitCompatibility()(*ih) ) { - result.emplace_back( std::move(ih) ); + auto est = estimator(&fl,es); + if (!est) return result; + using Algo = HitRCheck; + auto const & hitComp = (reinterpret_cast const&>(*est)).hitCompatibility(); + auto layerHits = layer.hits(); + result.reserve(layerHits.size()); + for (auto && ih : layerHits) { + if ( hitComp(*ih) ) + result.emplace_back( std::move(ih) ); } + } + } // std::cout << "RectangularEtaPhiTrackingRegion hits " << result.size() << std::endl; diff --git a/TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h b/TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h index b98c2724b90dc..7f56e595fa8fe 100755 --- a/TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h +++ b/TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h @@ -15,7 +15,7 @@ class EtaPhiMeasurementEstimator final : public MeasurementEstimator { public: - explicit EtaPhiMeasurementEstimator(double dEta, double dPhi) : + explicit EtaPhiMeasurementEstimator(float dEta, float dPhi) : thedEta(dEta), thedPhi(dPhi) {} @@ -34,8 +34,8 @@ class EtaPhiMeasurementEstimator final : public MeasurementEstimator { return new EtaPhiMeasurementEstimator(*this); } private: - double thedEta; - double thedPhi; + float thedEta; + float thedPhi; }; diff --git a/TrackingTools/KalmanUpdators/src/EtaPhiMeasurementEstimator.cc b/TrackingTools/KalmanUpdators/src/EtaPhiMeasurementEstimator.cc index f62adda36f7ce..468cc483d3ab2 100755 --- a/TrackingTools/KalmanUpdators/src/EtaPhiMeasurementEstimator.cc +++ b/TrackingTools/KalmanUpdators/src/EtaPhiMeasurementEstimator.cc @@ -10,34 +10,41 @@ std::pair EtaPhiMeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& aRecHit) const { - double dEta = fabs(tsos.globalPosition().eta() - aRecHit.globalPosition().eta()); - double dPhi = deltaPhi< double > (tsos.globalPosition().phi(), aRecHit.globalPosition().phi()); + // As the center of the plane is in the region, no need to waste time to test the hit as well + return std::make_pair(true, 1.0); + + /* HERE AS THE SECOND LINE WAS BUGGY... (AND MAY BE NEEDED IF WE MODIFY THE ESTIMATE WRT PLANE TO ACCOUNT FOR DET SPAN + + auto dEta = std::abs(tsos.globalPosition().eta() - aRecHit.globalPosition().eta()); + auto dPhi = std::abs(deltaPhi(tsos.globalPosition().barePhi(), aRecHit.globalPosition().barePhi())); LogDebug("EtaPhiMeasurementEstimator")<< " The state to compare with is \n"<< tsos << " The hit position is:\n" << aRecHit.globalPosition() << " deta: "<< dEta<< " dPhi: "< (tsos.globalPosition().phi(), plane.position().phi()); + auto dEta = std::abs(tsos.globalPosition().eta() - plane.eta()); + auto dPhi = std::abs(deltaPhi(tsos.globalPosition().barePhi(), plane.phi())); + LogDebug("EtaPhiMeasurementEstimator")<< "The state to compare with is \n"<< tsos << "\n" << "The plane position center is: " << plane.position() << "\n" << "the deta = " << thedEta << " --- the dPhi = " << thedPhi << "\n" << "deta = "<< fabs(dEta)<< " --- dPhi = "<