diff --git a/CommonTools/RecoAlgos/src/ClusterStorer.cc b/CommonTools/RecoAlgos/src/ClusterStorer.cc index e08c16c59fc07..760be3efc6c0b 100644 --- a/CommonTools/RecoAlgos/src/ClusterStorer.cc +++ b/CommonTools/RecoAlgos/src/ClusterStorer.cc @@ -9,6 +9,7 @@ #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" // FastSim hits: #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h" #include "DataFormats/TrackerRecHit2D/interface/FastProjectedTrackerRecHit.h" @@ -44,9 +45,12 @@ namespace helper { } else if (hit_type == typeid(Phase2TrackerRecHit1D)) { //FIXME:: this is just temporary solution for phase2, //it is not really running in the phase2 tracking wf - yet... - //std::cout << "| It is a Phase2TrackerRecHit1D hit !!" << std::endl; phase2OTClusterRecords_.push_back( Phase2OTClusterHitRecord(static_cast(newHit), hits, index)); + } else if (hit_type == typeid(VectorHit)) { + //FIXME:: this is just temporary solution for phase2, + //the VectorHit has 2 clusters but just a hit! + phase2OTClusterRecords_.push_back(Phase2OTClusterHitRecord(static_cast(newHit), hits, index)); } else { if (hit_type == typeid(FastTrackerRecHit) || hit_type == typeid(FastProjectedTrackerRecHit) || hit_type == typeid(FastMatchedTrackerRecHit)) { diff --git a/Configuration/ProcessModifiers/python/vectorHits_cff.py b/Configuration/ProcessModifiers/python/vectorHits_cff.py new file mode 100644 index 0000000000000..13c607fc1c922 --- /dev/null +++ b/Configuration/ProcessModifiers/python/vectorHits_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier is for activating Vector Hits reconstruction in the Phase II OT + +vectorHits = cms.Modifier() diff --git a/Configuration/PyReleaseValidation/python/relval_2026.py b/Configuration/PyReleaseValidation/python/relval_2026.py index e023aaf2653d8..3d753c1ea6c79 100644 --- a/Configuration/PyReleaseValidation/python/relval_2026.py +++ b/Configuration/PyReleaseValidation/python/relval_2026.py @@ -19,6 +19,7 @@ numWFIB.extend([23434.99,23434.999]) #2026D49 premixing combined stage1+stage2 (ttbar+PU200, ttbar+PU50 for PR test) numWFIB.extend([23234.21,23434.21]) #2026D49 prodlike, prodlike PU numWFIB.extend([23234.103]) #2026D49 aging +numWFIB.extend([23234.9]) #2026D49 vector hits numWFIB.extend([28234.0]) #2026D60 numWFIB.extend([29834.0]) #2026D64 numWFIB.extend([30234.0]) #2026D65 diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 0abe3e53bad20..0051cb40913b2 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -300,6 +300,23 @@ def condition_(self, fragment, stepList, key, hasHarvest): '--procModifiers': 'trackingMkFit' } +# Vector Hits workflows +class UpgradeWorkflow_vectorHits(UpgradeWorkflow): + def setup_(self, step, stepName, stepDict, k, properties): + stepDict[stepName][k] = merge([{'--procModifiers': 'vectorHits'}, stepDict[step][k]]) + def condition(self, fragment, stepList, key, hasHarvest): + return fragment=="TTbar_14TeV" and '2026' in key +upgradeWFs['vectorHits'] = UpgradeWorkflow_vectorHits( + steps = [ + 'RecoGlobal', + ], + PU = [ + 'RecoGlobal', + ], + suffix = '_vectorHits', + offset = 0.9, +) + # Patatrack workflows class UpgradeWorkflowPatatrack(UpgradeWorkflow): def condition(self, fragment, stepList, key, hasHarvest): diff --git a/DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h b/DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h index 9928e44849e48..3b55cce07b2ea 100644 --- a/DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h +++ b/DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h @@ -23,6 +23,8 @@ class BaseTrackerRecHit : public TrackingRecHit { // no position (as in persistent) BaseTrackerRecHit(DetId id, trackerHitRTTI::RTTI rt) : TrackingRecHit(id, (unsigned int)(rt)), qualWord_(0) {} + BaseTrackerRecHit(const GeomDet& idet, trackerHitRTTI::RTTI rt) + : TrackingRecHit(idet, (unsigned int)(rt)), qualWord_(0) {} BaseTrackerRecHit(const LocalPoint& p, const LocalError& e, GeomDet const& idet, trackerHitRTTI::RTTI rt) : TrackingRecHit(idet, (unsigned int)(rt)), pos_(p), err_(e), qualWord_(0) { @@ -49,18 +51,17 @@ class BaseTrackerRecHit : public TrackingRecHit { // verify that hits can share clusters... inline bool sameDetModule(TrackingRecHit const& hit) const; - bool hasPositionAndError() const final; + bool hasPositionAndError() const override; - LocalPoint localPosition() const final { + LocalPoint localPosition() const override { check(); return pos_; } - LocalError localPositionError() const final { + LocalError localPositionError() const override { check(); return err_; } - const LocalPoint& localPositionFast() const { check(); return pos_; diff --git a/DataFormats/TrackerRecHit2D/interface/TkCloner.h b/DataFormats/TrackerRecHit2D/interface/TkCloner.h index ed6ea47cd5d8d..cc98baef4ae8c 100644 --- a/DataFormats/TrackerRecHit2D/interface/TkCloner.h +++ b/DataFormats/TrackerRecHit2D/interface/TkCloner.h @@ -11,6 +11,7 @@ class SiStripRecHit1D; class SiStripMatchedRecHit2D; class ProjectedSiStripRecHit2D; class Phase2TrackerRecHit1D; +class VectorHit; class TkCloner { public: @@ -39,6 +40,7 @@ class TkCloner { TrajectoryStateOnSurface const& tsos) const = 0; virtual std::unique_ptr operator()(Phase2TrackerRecHit1D const& hit, TrajectoryStateOnSurface const& tsos) const = 0; + virtual std::unique_ptr operator()(VectorHit const& hit, TrajectoryStateOnSurface const& tsos) const = 0; #ifndef __GCCXML__ virtual TrackingRecHit::ConstRecHitPointer makeShared(SiPixelRecHit const& hit, @@ -53,6 +55,8 @@ class TkCloner { TrajectoryStateOnSurface const& tsos) const = 0; virtual TrackingRecHit::ConstRecHitPointer makeShared(Phase2TrackerRecHit1D const& hit, TrajectoryStateOnSurface const& tsos) const = 0; + virtual TrackingRecHit::ConstRecHitPointer makeShared(VectorHit const& hit, + TrajectoryStateOnSurface const& tsos) const = 0; #endif }; #endif diff --git a/DataFormats/TrackerRecHit2D/interface/VectorHit.h b/DataFormats/TrackerRecHit2D/interface/VectorHit.h new file mode 100644 index 0000000000000..df0a771f686d8 --- /dev/null +++ b/DataFormats/TrackerRecHit2D/interface/VectorHit.h @@ -0,0 +1,153 @@ +#ifndef DataFormats_TrackerRecHit2D_VectorHit_h +#define DataFormats_TrackerRecHit2D_VectorHit_h + +/** \class VectorHit + * + * 4-parameter RecHits for Phase2 Tracker (x,y, dx/dz, dy/dz) + * + * $Date: 2015/03/30 $ + * \author Erica Brondolin + * + */ + +#include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h" + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" + +#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "DataFormats/GeometryVector/interface/LocalVector.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h" + +#include "DataFormats/TrackerRecHit2D/interface/TkCloner.h" + +class VectorHit final : public BaseTrackerRecHit { +public: + typedef OmniClusterRef::Phase2Cluster1DRef ClusterRef; + + VectorHit() : thePosition(), theDirection(), theCovMatrix() { setType(bad); } + + VectorHit(const GeomDet& idet, + const LocalPoint& posInner, + const LocalVector& dir, + const AlgebraicSymMatrix44& covMatrix, + const float chi2, + OmniClusterRef const& lower, + OmniClusterRef const& upper, + const float curvature, + const float curvatureError, + const float phi); + + VectorHit(const GeomDet& idet, + const VectorHit2D& vh2Dzx, + const VectorHit2D& vh2Dzy, + OmniClusterRef const& lower, + OmniClusterRef const& upper, + const float curvature, + const float curvatureError, + const float phi); + + ~VectorHit() override = default; + + VectorHit* clone() const override { return new VectorHit(*this); } + RecHitPointer cloneSH() const override { return std::make_shared(*this); } + + bool sharesInput(const TrackingRecHit* other, SharedInputType what) const override; + bool sharesClusters(VectorHit const& other, SharedInputType what) const; + + // Parameters of the segment, for the track fit + // For a 4D segment: (dx/dz,dy/dz,x,y) + bool hasPositionAndError() const override { + //if det is present pos&err are available as well. + //if det() is not present (null) the hit has been read from file and not updated + return det(); + }; + + void getKfComponents(KfComponentsHolder& holder) const override { getKfComponents4D(holder); } + void getKfComponents4D(KfComponentsHolder& holder) const; + + // returning methods + LocalPoint localPosition() const override { return thePosition; } + virtual LocalVector localDirection() const { return theDirection; } + const AlgebraicSymMatrix44& covMatrix() const; + LocalError localPositionError() const override; + LocalError localDirectionError() const; + Global3DVector globalDirectionVH() const; + + float chi2() const { return theChi2; } + int dimension() const override { return theDimension; } + float curvature() const { return theCurvature; } + float curvatureError() const { return theCurvatureError; } + float phi() const { return thePhi; } + + float transverseMomentum(float magField) const; + float momentum(float magField) const; + + /// "lower" is logical, not geometrically lower; in pixel-strip modules the "lower" is always a pixel + ClusterRef lowerCluster() const { return theLowerCluster.cluster_phase2OT(); } + ClusterRef upperCluster() const { return theUpperCluster.cluster_phase2OT(); } + OmniClusterRef const lowerClusterRef() const { return theLowerCluster; } + OmniClusterRef const upperClusterRef() const { return theUpperCluster; } + + //FIXME::to update with a proper CPE maybe... + Global3DPoint lowerGlobalPos() const; + Global3DPoint upperGlobalPos() const; + static Global3DPoint phase2clusterGlobalPos(const PixelGeomDetUnit* geomDet, ClusterRef cluster); + GlobalError lowerGlobalPosErr() const; + GlobalError upperGlobalPosErr() const; + static GlobalError phase2clusterGlobalPosErr(const PixelGeomDetUnit* geomDet); + + bool isPhase2() const override { return true; } + + //FIXME: I have always two clusters in a VH + OmniClusterRef const& firstClusterRef() const override { return theLowerCluster; } + ClusterRef cluster() const { return theLowerCluster.cluster_phase2OT(); } + + //This method returns the direction of the segment/stub in global coordinates + Global3DVector globalDirection() const; + float theta() const; + + // Access to component RecHits (if any) + std::vector recHits() const override; + std::vector recHits() override; + +private: + // double dispatch + VectorHit* clone_(TkCloner const& cloner, TrajectoryStateOnSurface const& tsos) const override { + return cloner(*this, tsos).release(); + } + RecHitPointer cloneSH_(TkCloner const& cloner, TrajectoryStateOnSurface const& tsos) const override { + return cloner.makeShared(*this, tsos); + } + + LocalPoint thePosition; + LocalVector theDirection; + + // the covariance matrix, has the following meaning + // mat[0][0]=var(dx/dz) + // mat[1][1]=var(dy/dz) + // mat[2][2]=var(x) + // mat[3][3]=var(y) + // mat[0][2]=cov(dx/dz,x) + // mat[1][3]=cov(dy/dz,y) + AlgebraicSymMatrix44 theCovMatrix; + float theChi2; + static constexpr int theDimension = 4; + OmniClusterRef theLowerCluster; + OmniClusterRef theUpperCluster; + float theCurvature; + float theCurvatureError; + float thePhi; +}; + +inline bool operator<(const VectorHit& one, const VectorHit& other) { return (one.chi2() < other.chi2()); } + +std::ostream& operator<<(std::ostream& os, const VectorHit& vh); + +typedef edmNew::DetSetVector VectorHitCollection; + +#endif diff --git a/DataFormats/TrackerRecHit2D/interface/VectorHit2D.h b/DataFormats/TrackerRecHit2D/interface/VectorHit2D.h new file mode 100644 index 0000000000000..d55130f6df8f3 --- /dev/null +++ b/DataFormats/TrackerRecHit2D/interface/VectorHit2D.h @@ -0,0 +1,36 @@ +#ifndef DataFormats_TrackerRecHit2D_VectorHit2D_h +#define DataFormats_TrackerRecHit2D_VectorHit2D_h + +#include "DataFormats/GeometryVector/interface/LocalVector.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/CLHEP/interface/AlgebraicObjects.h" +#include "DataFormats/GeometrySurface/interface/LocalError.h" + +//Stores geometric information about VectorHits, their convariance matrix and chi2 of the compatability of the two hits + +class VectorHit2D { +public: + VectorHit2D() : thePosition(), theDirection(), theCovMatrix(), theChi2() {} + VectorHit2D(const LocalPoint& pos, const LocalVector& dir, const AlgebraicSymMatrix22& covMatrix, const double& chi2) + : thePosition(pos), + theDirection(dir), + theCovMatrix(covMatrix), + theLocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]), + theChi2(chi2) {} + + const LocalPoint& localPosition() const { return thePosition; } + const LocalVector& localDirection() const { return theDirection; } + const LocalError& localDirectionError() const { return theLocalError; } + const AlgebraicSymMatrix22& covMatrix() const { return theCovMatrix; } + float chi2() const { return theChi2; } + int dimension() const { return theDimension; } + +private: + LocalPoint thePosition; + LocalVector theDirection; + AlgebraicSymMatrix22 theCovMatrix; + LocalError theLocalError; + float theChi2; + static constexpr int theDimension = 2; +}; +#endif diff --git a/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h b/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h index 23370e4284049..b26198b498f69 100644 --- a/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h +++ b/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h @@ -17,7 +17,8 @@ namespace trackerHitRTTI { fastProjMono = 8, fastMatch = 9, notFromCluster = 10, - mipTiming = 11 + mipTiming = 11, + vector = 12 }; inline RTTI rtti(TrackingRecHit const& hit) { return RTTI(hit.getRTTI()); } inline bool isUndef(TrackingRecHit const& hit) { return rtti(hit) == undef; } @@ -32,10 +33,13 @@ namespace trackerHitRTTI { inline bool isMatched(TrackingRecHit const& hit) { return rtti(hit) == match || rtti(hit) == fastMatch; } inline bool isMulti(TrackingRecHit const& hit) { return rtti(hit) == multi; } inline bool isSingleType(TrackingRecHit const& hit) { return (rtti(hit) > 0) & (rtti(hit) < 4); } - inline bool isFromDet(TrackingRecHit const& hit) { return (rtti(hit) > 0) & (rtti(hit) < 6); } + inline bool isFromDet(TrackingRecHit const& hit) { return (((rtti(hit) > 0) & (rtti(hit) < 6)) | (rtti(hit) == 12)); } inline bool isFast(TrackingRecHit const& hit) { return (rtti(hit) > 5) & (rtti(hit) <= 9); } - inline bool isFromDetOrFast(TrackingRecHit const& hit) { return (rtti(hit) > 0) & (rtti(hit) < 10); } + inline bool isFromDetOrFast(TrackingRecHit const& hit) { + return (((rtti(hit) > 0) & (rtti(hit) < 10)) | (rtti(hit) == 12)); + } inline bool isTiming(TrackingRecHit const& hit) { return rtti(hit) == mipTiming; } + inline bool isVector(TrackingRecHit const& hit) { return rtti(hit) == vector; } inline unsigned int projId(TrackingRecHit const& hit) { return hit.rawId() + int(rtti(hit)) - 1; } } // namespace trackerHitRTTI diff --git a/DataFormats/TrackerRecHit2D/src/BaseTrackerRecHit.cc b/DataFormats/TrackerRecHit2D/src/BaseTrackerRecHit.cc index 4b1d7c60c8051..35aca8a3339c8 100644 --- a/DataFormats/TrackerRecHit2D/src/BaseTrackerRecHit.cc +++ b/DataFormats/TrackerRecHit2D/src/BaseTrackerRecHit.cc @@ -27,10 +27,9 @@ void BaseTrackerRecHit::check() const { #endif bool BaseTrackerRecHit::hasPositionAndError() const { + //if det is present pos&err are available as well. + // //if det() is not present (null) the hit has been read from file and not updated return det(); - - // return (err_.xx() != 0) || (err_.yy() != 0) || (err_.xy() != 0) || - // (pos_.x() != 0) || (pos_.y() != 0) || (pos_.z() != 0); } void BaseTrackerRecHit::getKfComponents1D(KfComponentsHolder &holder) const { diff --git a/DataFormats/TrackerRecHit2D/src/VectorHit.cc b/DataFormats/TrackerRecHit2D/src/VectorHit.cc new file mode 100644 index 0000000000000..d57358cc8f5d0 --- /dev/null +++ b/DataFormats/TrackerRecHit2D/src/VectorHit.cc @@ -0,0 +1,188 @@ +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "Geometry/CommonDetUnit/interface/StackGeomDet.h" +#include "CLHEP/Units/PhysicalConstants.h" + +VectorHit::VectorHit(const GeomDet& idet, + const LocalPoint& posLower, + const LocalVector& dir, + const AlgebraicSymMatrix44& covMatrix, + const float chi2, + OmniClusterRef const& lower, + OmniClusterRef const& upper, + const float curvature, + const float curvatureError, + const float phi) + : BaseTrackerRecHit(idet, trackerHitRTTI::vector), + thePosition(posLower), + theDirection(dir), + theCovMatrix(covMatrix), + theChi2(chi2), + theLowerCluster(lower), + theUpperCluster(upper), + theCurvature(curvature), + theCurvatureError(curvatureError), + thePhi(phi) {} + +VectorHit::VectorHit(const GeomDet& idet, + const VectorHit2D& vh2Dzx, + const VectorHit2D& vh2Dzy, + OmniClusterRef const& lower, + OmniClusterRef const& upper, + const float curvature, + const float curvatureError, + const float phi) + : BaseTrackerRecHit(idet, trackerHitRTTI::vector), + thePosition(LocalPoint(vh2Dzx.localPosition().x(), vh2Dzy.localPosition().x(), 0.)), + theDirection(LocalVector(vh2Dzx.localDirection().x(), vh2Dzy.localDirection().x(), 1.)), + theLowerCluster(lower), + theUpperCluster(upper), + theCurvature(curvature), + theCurvatureError(curvatureError), + thePhi(phi) { + //building the cov matrix 4x4 starting from the 2x2 + const AlgebraicSymMatrix22& covMatZX = vh2Dzx.covMatrix(); + const AlgebraicSymMatrix22& covMatZY = vh2Dzy.covMatrix(); + + theCovMatrix = AlgebraicSymMatrix44(); + + theCovMatrix[0][0] = covMatZX[0][0]; // var(dx/dz) + theCovMatrix[1][1] = covMatZY[0][0]; // var(dy/dz) + theCovMatrix[2][2] = covMatZX[1][1]; // var(x) + theCovMatrix[3][3] = covMatZY[1][1]; // var(y) + theCovMatrix[0][2] = covMatZX[0][1]; // cov(dx/dz,x) + theCovMatrix[1][3] = covMatZY[0][1]; // cov(dy/dz,y) + + theChi2 = vh2Dzx.chi2() + vh2Dzy.chi2(); +} + +bool VectorHit::sharesInput(const TrackingRecHit* other, SharedInputType what) const { + if (what == all && (geographicalId() != other->geographicalId())) + return false; + + if (!sameDetModule(*other)) + return false; + + if (trackerHitRTTI::isVector(*other)) { + const VectorHit* otherVh = static_cast(other); + return sharesClusters(*otherVh, what); + } + + if (what == all) + return false; + + // what about multi??? + auto const& otherClus = reinterpret_cast(other)->firstClusterRef(); + return (otherClus == lowerClusterRef()) || (otherClus == upperClusterRef()); +} + +bool VectorHit::sharesClusters(VectorHit const& other, SharedInputType what) const { + bool lower = this->lowerClusterRef() == other.lowerClusterRef(); + bool upper = this->upperClusterRef() == other.upperClusterRef(); + + return (what == TrackingRecHit::all) ? (lower && upper) : (upper || lower); +} + +void VectorHit::getKfComponents4D(KfComponentsHolder& holder) const { + AlgebraicVector4& pars = holder.params(); + pars[0] = theDirection.x(); + pars[1] = theDirection.y(); + pars[2] = thePosition.x(); + pars[3] = thePosition.y(); + + holder.errors() = theCovMatrix; + + ProjectMatrix& pf = holder.projFunc(); + for (int i = 0; i < 4; ++i) + pf.index[i] = i + 1; + + holder.measuredParams() = AlgebraicVector4(&holder.tsosLocalParameters().At(1), theDimension); + holder.measuredErrors() = holder.tsosLocalErrors().Sub(1, 1); +} + +Global3DPoint VectorHit::lowerGlobalPos() const { + const StackGeomDet* stackDet = static_cast(det()); + const PixelGeomDetUnit* geomDetLower = static_cast(stackDet->lowerDet()); + return phase2clusterGlobalPos(geomDetLower, lowerCluster()); +} + +Global3DPoint VectorHit::upperGlobalPos() const { + const StackGeomDet* stackDet = static_cast(det()); + const PixelGeomDetUnit* geomDetUpper = static_cast(stackDet->upperDet()); + return phase2clusterGlobalPos(geomDetUpper, upperCluster()); +} + +Global3DPoint VectorHit::phase2clusterGlobalPos(const PixelGeomDetUnit* geomDet, ClusterRef cluster) { + const PixelTopology* topo = &geomDet->specificTopology(); + float ix = cluster->center(); + float iy = cluster->column() + 0.5f; // halfway the column + LocalPoint lp(topo->localX(ix), topo->localY(iy), 0); // x, y, z + Global3DPoint gp = geomDet->surface().toGlobal(lp); + return gp; +} + +GlobalError VectorHit::lowerGlobalPosErr() const { + const StackGeomDet* stackDet = static_cast(det()); + const PixelGeomDetUnit* geomDetLower = static_cast(stackDet->lowerDet()); + return phase2clusterGlobalPosErr(geomDetLower); +} + +GlobalError VectorHit::upperGlobalPosErr() const { + const StackGeomDet* stackDet = static_cast(det()); + const PixelGeomDetUnit* geomDetUpper = static_cast(stackDet->upperDet()); + return phase2clusterGlobalPosErr(geomDetUpper); +} + +GlobalError VectorHit::phase2clusterGlobalPosErr(const PixelGeomDetUnit* geomDet) { + const PixelTopology* topo = &geomDet->specificTopology(); + float pitchX = topo->pitch().first; + float pitchY = topo->pitch().second; + constexpr float invTwelve = 1. / 12; + LocalError le(pow(pitchX, 2) * invTwelve, 0, pow(pitchY, 2) * invTwelve); // e2_xx, e2_xy, e2_yy + GlobalError ge(ErrorFrameTransformer().transform(le, geomDet->surface())); + return ge; +} + +Global3DVector VectorHit::globalDirectionVH() const { + Local3DVector theLocalDelta = + LocalVector(theDirection.x() * theDirection.z(), theDirection.y() * theDirection.z(), theDirection.z()); + Global3DVector g = det()->surface().toGlobal(theLocalDelta); + return g; +} + +Global3DVector VectorHit::globalDirection() const { return (det()->surface().toGlobal(localDirection())); } + +float VectorHit::theta() const { return globalDirection().theta(); } + +float VectorHit::transverseMomentum(float magField) const { + return magField * (CLHEP::c_light * 1e-5F) / theCurvature; + // pT [GeV] ~ 0.3 * B[T] * R [m], curvature is in cms, using precise value from speed of light + // because curvature is a signed quantity this transverse momentum is also signed +} +float VectorHit::momentum(float magField) const { return transverseMomentum(magField) / (1. * sin(theta())); } + +LocalError VectorHit::localPositionError() const { + return LocalError(theCovMatrix[2][2], theCovMatrix[2][3], theCovMatrix[3][3]); +} + +LocalError VectorHit::localDirectionError() const { + return LocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]); +} + +const AlgebraicSymMatrix44& VectorHit::covMatrix() const { return theCovMatrix; } + +std::ostream& operator<<(std::ostream& os, const VectorHit& vh) { + os << " VectorHit create in the DetId#: " << vh.geographicalId() << "\n" + << " Vectorhit local position : " << vh.localPosition() << "\n" + << " Vectorhit local direction : " << vh.localDirection() << "\n" + << " Vectorhit global direction : " << vh.globalDirection() << "\n" + << " Lower cluster global position : " << vh.lowerGlobalPos() << "\n" + << " Upper cluster global position : " << vh.upperGlobalPos(); + + return os; +} + +/// Access to component RecHits (if any) +std::vector VectorHit::recHits() const { return {}; } + +/// Non-const access to component RecHits (if any) +std::vector VectorHit::recHits() { return {}; } diff --git a/DataFormats/TrackerRecHit2D/src/classes.h b/DataFormats/TrackerRecHit2D/src/classes.h index 45302da925815..03140aa37a5bb 100644 --- a/DataFormats/TrackerRecHit2D/src/classes.h +++ b/DataFormats/TrackerRecHit2D/src/classes.h @@ -12,6 +12,7 @@ #include "DataFormats/Common/interface/RefProd.h" #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" #include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" #include "DataFormats/Common/interface/Wrapper.h" #include "DataFormats/Common/interface/Ref.h" #include "DataFormats/Common/interface/RefVector.h" @@ -27,6 +28,7 @@ #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" #include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include #endif // SISTRIPRECHIT_CLASSES_H diff --git a/DataFormats/TrackerRecHit2D/src/classes_def.xml b/DataFormats/TrackerRecHit2D/src/classes_def.xml index e2a6415e917ad..5de054ebdfa8f 100644 --- a/DataFormats/TrackerRecHit2D/src/classes_def.xml +++ b/DataFormats/TrackerRecHit2D/src/classes_def.xml @@ -189,5 +189,18 @@ + + + + + + + + + + + + + diff --git a/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc b/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc index a58bc1841e0af..55dd2b5f9458e 100644 --- a/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc +++ b/Geometry/TrackerGeometryBuilder/src/TrackerGeomBuilderFromGeometricDet.cc @@ -262,10 +262,10 @@ void TrackerGeomBuilderFromGeometricDet::buildGeomDet(TrackerGeometry* tracker) tracker->addDetId(composedDetId); } else if (gduTypeName.find("Lower") != std::string::npos) { - //FIXME::ERICA: the plane builder is built in the middle... - PlaneBuilderForGluedDet::ResultType plane = gluedplaneBuilder.plane(composed); + //The plane is *not* built in the middle, but on the Lower surface + Plane* plane = new Plane(dus->surface()); composedDetId = theTopo->stack(gduId[i]); - StackGeomDet* stackDet = new StackGeomDet(&(*plane), dum, dus, composedDetId); + StackGeomDet* stackDet = new StackGeomDet(&(*plane), dus, dum, composedDetId); tracker->addDet((GeomDet*)stackDet); tracker->addDetId(composedDetId); } diff --git a/RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py b/RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py index 4d25967bd803f..3cae176059b3b 100644 --- a/RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py +++ b/RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py @@ -23,6 +23,7 @@ from RecoLocalTracker.SiPhase2Clusterizer.phase2TrackerClusterizer_cfi import * from RecoLocalTracker.Phase2TrackerRecHits.Phase2StripCPEGeometricESProducer_cfi import * +from RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2RecHitMatcher_cfi import * _pixeltrackerlocalrecoTask_phase2 = pixeltrackerlocalrecoTask.copy() _pixeltrackerlocalrecoTask_phase2.add(siPhase2Clusters) diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/BuildFile.xml b/RecoLocalTracker/SiPhase2VectorHitBuilder/BuildFile.xml new file mode 100644 index 0000000000000..67dc8d580bb3d --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/BuildFile.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h b/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h new file mode 100644 index 0000000000000..8d31827293481 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h @@ -0,0 +1,89 @@ +//--------------------------------------------------------------------------- +// class VectorHitBuilderAlgorithm +// author: ebrondol,nathera +// date: May, 2015 +//--------------------------------------------------------------------------- + +#ifndef RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitBuilderAlgorithm_H +#define RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitBuilderAlgorithm_H + +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "CommonTools/Statistics/interface/LinearFit.h" + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +class VectorHitBuilderAlgorithm : public VectorHitBuilderAlgorithmBase { +public: + VectorHitBuilderAlgorithm(const edm::ParameterSet& conf, + const TrackerGeometry* tkGeomProd, + const TrackerTopology* tkTopoProd, + const ClusterParameterEstimator* cpeProd) + : VectorHitBuilderAlgorithmBase(conf, tkGeomProd, tkTopoProd, cpeProd){}; + ~VectorHitBuilderAlgorithm() override = default; + + void run(edm::Handle> clusters, + VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + edmNew::DetSetVector& clustersAcc, + edmNew::DetSetVector& clustersRej) const override; + + //not implemented yet + bool checkClustersCompatibilityBeforeBuilding(edm::Handle> clusters, + const Detset& theLowerDetSet, + const Detset& theUpperDetSet) const; + bool checkClustersCompatibility(Local3DPoint& posinner, + Local3DPoint& posouter, + LocalError& errinner, + LocalError& errouter) const; + struct CurvatureAndPhi { + float curvature; + float curvatureError; + float phi; + }; + CurvatureAndPhi curvatureANDphi(Global3DPoint gPositionLower, + Global3DPoint gPositionUpper, + GlobalError gErrorLower, + GlobalError gErrorUpper) const; + + void buildVectorHits(VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + DetId detIdStack, + const StackGeomDet* stack, + edm::Handle> clusters, + const Detset& DSVinner, + const Detset& DSVouter, + const std::vector& phase2OTClustersToSkip = std::vector()) const override; + + VectorHit buildVectorHit(const StackGeomDet* stack, + Phase2TrackerCluster1DRef lower, + Phase2TrackerCluster1DRef upper) const override; + + void fit2Dzx(const Local3DPoint lpCI, + const Local3DPoint lpCO, + const LocalError leCI, + const LocalError leCO, + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const; + void fit2Dzy(const Local3DPoint lpCI, + const Local3DPoint lpCO, + const LocalError leCI, + const LocalError leCO, + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const; + + void fit(float x[2], + float y[2], + float sigy[2], + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const; +}; + +#endif diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h b/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h new file mode 100644 index 0000000000000..61d80a95cff87 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h @@ -0,0 +1,71 @@ +#ifndef RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitBuilderAlgorithmBase_H +#define RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitBuilderAlgorithmBase_H + +#include "FWCore/Framework/interface/ESHandle.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonDetUnit/interface/StackGeomDet.h" +#include "RecoLocalTracker/Phase2TrackerRecHits/interface/Phase2StripCPE.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" +#include "FWCore/Framework/interface/FrameworkfwdMostUsed.h" + +class VectorHitBuilderAlgorithmBase { +public: + typedef edm::Ref, Phase2TrackerCluster1D> Phase2TrackerCluster1DRef; + typedef edmNew::DetSet Detset; + typedef Detset::const_iterator const_iterator; + typedef edmNew::DetSetVector output_t; + typedef std::pair> StackClusters; + + VectorHitBuilderAlgorithmBase(const edm::ParameterSet&, + const TrackerGeometry*, + const TrackerTopology*, + const ClusterParameterEstimator*); + virtual ~VectorHitBuilderAlgorithmBase() {} + + //FIXME::ERICA::this should be template, return different collection for different algo used!! + virtual void run(edm::Handle> clusters, + VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + edmNew::DetSetVector& clustersAcc, + edmNew::DetSetVector& clustersRej) const = 0; + + virtual void buildVectorHits(VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + DetId detIdStack, + const StackGeomDet* stack, + edm::Handle> clusters, + const Detset& DSVinner, + const Detset& DSVouter, + const std::vector& phase2OTClustersToSkip = std::vector()) const = 0; + + virtual VectorHit buildVectorHit(const StackGeomDet* stack, + Phase2TrackerCluster1DRef lower, + Phase2TrackerCluster1DRef upper) const = 0; + + double computeParallaxCorrection(const PixelGeomDetUnit*, + const Point3DBase&, + const PixelGeomDetUnit*, + const Point3DBase&) const; + + void printClusters(const edmNew::DetSetVector& clusters) const; + void printCluster(const GeomDet* geomDetUnit, const Phase2TrackerCluster1D* cluster) const; + + const TrackerGeometry* tkGeom_; + const TrackerTopology* tkTopo_; + const ClusterParameterEstimator* cpe_; + unsigned int nMaxVHforeachStack_; + std::vector barrelCut_; + std::vector endcapCut_; + +private: + edm::ESInputTag cpeTag_; +}; + +#endif diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/BuildFile.xml b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/BuildFile.xml new file mode 100644 index 0000000000000..1e43153bea1ba --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/SiPhase2RecHitMatcherESProducer.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/SiPhase2RecHitMatcherESProducer.cc new file mode 100644 index 0000000000000..974a59fc59682 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/SiPhase2RecHitMatcherESProducer.cc @@ -0,0 +1,84 @@ + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/Utilities/interface/EDMException.h" +#include "RecoLocalTracker/Records/interface/TkPhase2OTCPERecord.h" +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" + +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" +#include + +class SiPhase2RecHitMatcherESProducer : public edm::ESProducer { +public: + SiPhase2RecHitMatcherESProducer(const edm::ParameterSet&); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + std::unique_ptr produce(const TkPhase2OTCPERecord&); + +private: + std::string name_; + edm::ParameterSet pset_; + edm::ESGetToken geometryToken_; + edm::ESGetToken trackerTopoToken_; + edm::ESGetToken, TkPhase2OTCPERecord> cpeToken_; +}; + +SiPhase2RecHitMatcherESProducer::SiPhase2RecHitMatcherESProducer(const edm::ParameterSet& p) { + name_ = p.getParameter("ComponentName"); + if (!(name_ == "SiPhase2VectorHitMatcher")) { + throw cms::Exception("ConfigurationError") << "Configuration specifies unknown ComponentName .\n" + << "Currently only 'SiPhase2VectorHitMatcher' is supported\n"; + } + pset_ = p; + auto cc = setWhatProduced(this, name_); + geometryToken_ = cc.consumes(); + trackerTopoToken_ = cc.consumes(); + cpeToken_ = cc.consumes(p.getParameter("CPE")); +} + +std::unique_ptr SiPhase2RecHitMatcherESProducer::produce(const TkPhase2OTCPERecord& iRecord) { + std::unique_ptr matcher = std::make_unique( + pset_, + &iRecord.get(geometryToken_), + &iRecord.getRecord().get(trackerTopoToken_), + &iRecord.get(cpeToken_)); + + return matcher; +} + +void SiPhase2RecHitMatcherESProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("offlinestubs", "vectorHits"); + desc.add("maxVectorHits", 999999999); + desc.add("Algorithm", "VectorHitBuilderAlgorithm"); + desc.add("ComponentName", "SiPhase2VectorHitMatcher"); + desc.add("CPE", edm::ESInputTag("", "Phase2StripCPE")); + desc.add>("BarrelCut", + { + 0.0, + 0.05, + 0.06, + 0.08, + 0.09, + 0.12, + 0.2, + }); + desc.add("Phase2CPE_name", "Phase2StripCPE"); + desc.add("Clusters", "siPhase2Clusters"); + desc.add("maxVectorHitsInAStack", 999); + desc.add>("EndcapCut", + { + 0.0, + 0.1, + 0.1, + 0.1, + 0.1, + 0.1, + }); + descriptions.add("siPhase2RecHitMatcher", desc); +} + +DEFINE_FWK_EVENTSETUP_MODULE(SiPhase2RecHitMatcherESProducer); diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/VectorHitBuilderEDProducer.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/VectorHitBuilderEDProducer.cc new file mode 100644 index 0000000000000..404ce82f97004 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/plugins/VectorHitBuilderEDProducer.cc @@ -0,0 +1,125 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h" +#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "RecoLocalTracker/Records/interface/TkPhase2OTCPERecord.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +class VectorHitBuilderEDProducer : public edm::stream::EDProducer<> { +public: + explicit VectorHitBuilderEDProducer(const edm::ParameterSet&); + ~VectorHitBuilderEDProducer() override = default; + void produce(edm::Event&, const edm::EventSetup&) override; + void run(edm::Handle> clusters, + edmNew::DetSetVector& clustersAcc, + edmNew::DetSetVector& clustersRej, + VectorHitCollection& outputAcc, + VectorHitCollection& outputRej); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + const VectorHitBuilderAlgorithm* algo() const { return stubsBuilder_; } + +private: + const VectorHitBuilderAlgorithm* stubsBuilder_; + std::string offlinestubsTag_; + unsigned int maxOfflinestubs_; + edm::ESGetToken stubsBuilderToken_; + edm::EDGetTokenT> clusterProducer_; +}; + +VectorHitBuilderEDProducer::VectorHitBuilderEDProducer(edm::ParameterSet const& conf) + : offlinestubsTag_(conf.getParameter("offlinestubs")), + maxOfflinestubs_(conf.getParameter("maxVectorHits")), + stubsBuilderToken_(esConsumes(conf.getParameter("Algorithm"))) { + clusterProducer_ = + consumes>(edm::InputTag(conf.getParameter("Clusters"))); + + produces>("accepted"); + produces>("rejected"); + produces("accepted"); + produces("rejected"); +} + +void VectorHitBuilderEDProducer::produce(edm::Event& event, const edm::EventSetup& es) { + LogDebug("VectorHitBuilderEDProducer") << "VectorHitBuilderEDProducer::produce() begin"; + + // get input clusters data + auto clustersHandle = event.getHandle(clusterProducer_); + // create the final output collection + auto outputClustersAccepted = std::make_unique>(); + auto outputClustersRejected = std::make_unique>(); + std::unique_ptr outputVHAccepted(new VectorHitCollection()); + std::unique_ptr outputVHRejected(new VectorHitCollection()); + + stubsBuilder_ = &es.getData(stubsBuilderToken_); +#ifdef EDM_ML_DEBUG + // check on the input clusters + stubsBuilder_->printClusters(*clustersHandle); +#endif //EDM_ML_DEBUG + + // running the stub building algorithm + //ERICA::output should be moved in the different algo classes? + run(clustersHandle, *outputClustersAccepted, *outputClustersRejected, *outputVHAccepted, *outputVHRejected); +#ifdef EDM_ML_DEBUG + unsigned int numberOfVectorHits = 0; + for (const auto& dSViter : *outputVHAccepted) { + for (const auto& vh : dSViter) { + numberOfVectorHits++; + LogDebug("VectorHitBuilderEDProducer") << "\t vectorhit in output " << vh; + } + } + LogDebug("VectorHitBuilderEDProducer") << "found\n" << numberOfVectorHits << " .\n"; +#endif //EDM_ML_DEBUG + // write output to file + event.put(std::move(outputClustersAccepted), "accepted"); + event.put(std::move(outputClustersRejected), "rejected"); + event.put(std::move(outputVHAccepted), "accepted"); + event.put(std::move(outputVHRejected), "rejected"); +} + +void VectorHitBuilderEDProducer::run(edm::Handle> clusters, + edmNew::DetSetVector& clustersAcc, + edmNew::DetSetVector& clustersRej, + VectorHitCollection& outputAcc, + VectorHitCollection& outputRej) { + stubsBuilder_->run(clusters, outputAcc, outputRej, clustersAcc, clustersRej); +} +void VectorHitBuilderEDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("offlinestubs", "vectorHits"); + desc.add("maxVectorHits", 999999999); + desc.add("Algorithm", edm::ESInputTag("", "SiPhase2VectorHitMatcher")); + desc.add("CPE", edm::ESInputTag("phase2StripCPEESProducer", "Phase2StripCPE")); + desc.add>("BarrelCut", + { + 0.0, + 0.05, + 0.06, + 0.08, + 0.09, + 0.12, + 0.2, + }); + desc.add("Clusters", "siPhase2Clusters"); + desc.add("maxVectorHitsInAStack", 999); + desc.add>("EndcapCut", + { + 0.0, + 0.1, + 0.1, + 0.1, + 0.1, + 0.1, + }); + descriptions.add("siPhase2VectorHits", desc); +} +DEFINE_FWK_MODULE(VectorHitBuilderEDProducer); +#include "FWCore/Utilities/interface/typelookup.h" +TYPELOOKUP_DATA_REG(VectorHitBuilderEDProducer); diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/src/ES_VectorHitBuilderAlgorithm.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/ES_VectorHitBuilderAlgorithm.cc new file mode 100644 index 0000000000000..7bbb14a2b4f16 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/ES_VectorHitBuilderAlgorithm.cc @@ -0,0 +1,4 @@ +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" +#include "FWCore/Utilities/interface/typelookup.h" + +TYPELOOKUP_DATA_REG(VectorHitBuilderAlgorithm); diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithm.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithm.cc new file mode 100644 index 0000000000000..be41bfc3aa412 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithm.cc @@ -0,0 +1,514 @@ +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/ClusterParameterEstimator.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonTopologies/interface/PixelTopology.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit2D.h" + +void VectorHitBuilderAlgorithm::run(edm::Handle> clusters, + VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + edmNew::DetSetVector& clustersAcc, + edmNew::DetSetVector& clustersRej) const { + LogDebug("VectorHitBuilderAlgorithm") << "Run VectorHitBuilderAlgorithm ... \n"; + const auto* clustersPhase2Collection = clusters.product(); + + //loop over the DetSetVector + LogDebug("VectorHitBuilderAlgorithm") << "with #clusters : " << clustersPhase2Collection->size() << std::endl; + for (auto dSViter : *clustersPhase2Collection) { + unsigned int rawDetId1(dSViter.detId()); + DetId detId1(rawDetId1); + DetId lowerDetId, upperDetId; + if (tkTopo_->isLower(detId1)) { + lowerDetId = detId1; + upperDetId = tkTopo_->partnerDetId(detId1); + } else + continue; + + DetId detIdStack = tkTopo_->stack(detId1); + + //debug + LogDebug("VectorHitBuilderAlgorithm") << " DetId stack : " << detIdStack.rawId() << std::endl; + LogDebug("VectorHitBuilderAlgorithm") << " DetId lower set of clusters : " << lowerDetId.rawId(); + LogDebug("VectorHitBuilderAlgorithm") << " DetId upper set of clusters : " << upperDetId.rawId() << std::endl; + + const GeomDet* gd; + const StackGeomDet* stackDet; + const auto& it_detLower = dSViter; + const auto& it_detUpper = clustersPhase2Collection->find(upperDetId); + + if (it_detUpper != clustersPhase2Collection->end()) { + gd = tkGeom_->idToDet(detIdStack); + stackDet = dynamic_cast(gd); + buildVectorHits(vhAcc, vhRej, detIdStack, stackDet, clusters, it_detLower, *it_detUpper); + } + } + LogDebug("VectorHitBuilderAlgorithm") << "End run VectorHitBuilderAlgorithm ... \n"; +} + +bool VectorHitBuilderAlgorithm::checkClustersCompatibilityBeforeBuilding( + edm::Handle> clusters, + const Detset& theLowerDetSet, + const Detset& theUpperDetSet) const { + if (theLowerDetSet.size() == 1 && theUpperDetSet.size() == 1) + return true; + + //order lower clusters in u + std::vector lowerClusters; + lowerClusters.reserve(theLowerDetSet.size()); + if (theLowerDetSet.size() > 1) + LogDebug("VectorHitBuilderAlgorithm") << " more than 1 lower cluster! " << std::endl; + if (theUpperDetSet.size() > 1) + LogDebug("VectorHitBuilderAlgorithm") << " more than 1 upper cluster! " << std::endl; + for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) { + Phase2TrackerCluster1DRef clusterLower = edmNew::makeRefTo(clusters, cil); + lowerClusters.push_back(*clusterLower); + } + return true; +} + +bool VectorHitBuilderAlgorithm::checkClustersCompatibility(Local3DPoint& poslower, + Local3DPoint& posupper, + LocalError& errlower, + LocalError& errupper) const { + return true; +} + +//---------------------------------------------------------------------------- +//ERICA::in the DT code the global position is used to compute the alpha angle and put a cut on that. +void VectorHitBuilderAlgorithm::buildVectorHits(VectorHitCollection& vhAcc, + VectorHitCollection& vhRej, + DetId detIdStack, + const StackGeomDet* stack, + edm::Handle> clusters, + const Detset& theLowerDetSet, + const Detset& theUpperDetSet, + const std::vector& phase2OTClustersToSkip) const { + if (checkClustersCompatibilityBeforeBuilding(clusters, theLowerDetSet, theUpperDetSet)) { + LogDebug("VectorHitBuilderAlgorithm") << " compatible -> continue ... " << std::endl; + } else { + LogTrace("VectorHitBuilderAlgorithm") << " not compatible, going to the next cluster"; + } + + edmNew::DetSetVector::FastFiller vh_colAcc(vhAcc, detIdStack); + edmNew::DetSetVector::FastFiller vh_colRej(vhRej, detIdStack); + + unsigned int layerStack = tkTopo_->layer(stack->geographicalId()); + if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB) + LogDebug("VectorHitBuilderAlgorithm") << " \t is barrel. " << std::endl; + if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC) + LogDebug("VectorHitBuilderAlgorithm") << " \t is endcap. " << std::endl; + LogDebug("VectorHitBuilderAlgorithm") << " \t layer is : " << layerStack << std::endl; + + float cut = 0.0; + if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTB) + cut = barrelCut_.at(layerStack); + if (stack->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC) + cut = endcapCut_.at(layerStack); + LogDebug("VectorHitBuilderAlgorithm") << " \t the cut is:" << cut << std::endl; + + //only cache local parameters for upper cluster as we loop over lower clusters only once anyway + std::vector> localParamsUpper; + std::vector localGDUUpper; + + const PixelGeomDetUnit* gduUpp = dynamic_cast(stack->upperDet()); + for (auto const& clusterUpper : theUpperDetSet) { + localGDUUpper.push_back(gduUpp); + localParamsUpper.push_back(cpe_->localParameters(clusterUpper, *gduUpp)); + } + int upperIterator = 0; + const PixelGeomDetUnit* gduLow = dynamic_cast(stack->lowerDet()); + for (const_iterator cil = theLowerDetSet.begin(); cil != theLowerDetSet.end(); ++cil) { + LogDebug("VectorHitBuilderAlgorithm") << " lower clusters " << std::endl; + Phase2TrackerCluster1DRef cluL = edmNew::makeRefTo(clusters, cil); +#ifdef EDM_ML_DEBUG + printCluster(stack->lowerDet(), &*cluL); +#endif + auto&& lparamsLow = cpe_->localParameters(*cluL, *gduLow); + upperIterator = 0; + for (const_iterator ciu = theUpperDetSet.begin(); ciu != theUpperDetSet.end(); ++ciu) { + LogDebug("VectorHitBuilderAlgorithm") << "\t upper clusters " << std::endl; + Phase2TrackerCluster1DRef cluU = edmNew::makeRefTo(clusters, ciu); +#ifdef EDM_ML_DEBUG + printCluster(stack->upperDet(), &*cluU); +#endif + //applying the parallax correction + double pC = computeParallaxCorrection( + gduLow, lparamsLow.first, localGDUUpper[upperIterator], localParamsUpper[upperIterator].first); + LogDebug("VectorHitBuilderAlgorithm") << " \t parallax correction:" << pC << std::endl; + double lpos_upp_corr = 0.0; + double lpos_low_corr = 0.0; + auto const localUpperX = localParamsUpper[upperIterator].first.x(); + if (localUpperX > lparamsLow.first.x()) { + if (localUpperX > 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = localParamsUpper[upperIterator].first.x() - std::abs(pC); + } else if (localUpperX < 0) { + lpos_low_corr = lparamsLow.first.x() + std::abs(pC); + lpos_upp_corr = localUpperX; + } + } else if (localUpperX < lparamsLow.first.x()) { + if (localUpperX > 0) { + lpos_low_corr = lparamsLow.first.x() - std::abs(pC); + lpos_upp_corr = localUpperX; + } else if (localUpperX < 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = localUpperX + std::abs(pC); + } + } else { + if (localUpperX > 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = localUpperX - std::abs(pC); + } else if (localUpperX < 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = localUpperX + std::abs(pC); + } + } + + LogDebug("VectorHitBuilderAlgorithm") << " \t local pos upper corrected (x):" << lpos_upp_corr << std::endl; + LogDebug("VectorHitBuilderAlgorithm") << " \t local pos lower corrected (x):" << lpos_low_corr << std::endl; + + double width = lpos_low_corr - lpos_upp_corr; + LogDebug("VectorHitBuilderAlgorithm") << " \t width: " << width << std::endl; + + //old cut: indipendent from layer + //building my tolerance : 10*sigma + //double delta = 10.0 * sqrt(lparamsLow.second.xx() + localParamsUpper[upperIterator].second.xx()); + //LogDebug("VectorHitBuilderAlgorithm") << " \t delta: " << delta << std::endl; + //if( (lpos_upp_corr < lpos_low_corr + delta) && + // (lpos_upp_corr > lpos_low_corr - delta) ){ + //new cut: dependent on layers + if (std::abs(width) < cut) { + LogDebug("VectorHitBuilderAlgorithm") << " accepting VH! " << std::endl; + VectorHit vh = buildVectorHit(stack, cluL, cluU); + //protection: the VH can also be empty!! + if (vh.isValid()) { + vh_colAcc.push_back(vh); + } + + } else { + LogDebug("VectorHitBuilderAlgorithm") << " rejecting VH: " << std::endl; + //storing vh rejected for combinatiorial studies + VectorHit vh = buildVectorHit(stack, cluL, cluU); + if (vh.isValid()) { + vh_colRej.push_back(vh); + } + } + upperIterator += 1; + } + } +} + +VectorHit VectorHitBuilderAlgorithm::buildVectorHit(const StackGeomDet* stack, + Phase2TrackerCluster1DRef lower, + Phase2TrackerCluster1DRef upper) const { + LogTrace("VectorHitBuilderAlgorithm") << "Build VH with: "; +#ifdef EDM_ML_DEBUG + printCluster(stack->upperDet(), &*upper); +#endif + const PixelGeomDetUnit* geomDetLower = static_cast(stack->lowerDet()); + const PixelGeomDetUnit* geomDetUpper = static_cast(stack->upperDet()); + + auto&& lparamsLower = cpe_->localParameters(*lower, *geomDetLower); // x, y, z, e2_xx, e2_xy, e2_yy + Global3DPoint gparamsLower = geomDetLower->surface().toGlobal(lparamsLower.first); + LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower; + + auto&& lparamsUpper = cpe_->localParameters(*upper, *geomDetUpper); + Global3DPoint gparamsUpper = geomDetUpper->surface().toGlobal(lparamsUpper.first); + LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper; + + //local parameters of upper cluster in lower system of reference + Local3DPoint lparamsUpperInLower = geomDetLower->surface().toLocal(gparamsUpper); + + LogTrace("VectorHitBuilderAlgorithm") << "\t lower global pos: " << gparamsLower; + LogTrace("VectorHitBuilderAlgorithm") << "\t upper global pos: " << gparamsUpper; + + LogTrace("VectorHitBuilderAlgorithm") << "A:\t lower local pos: " << lparamsLower.first + << " with error: " << lparamsLower.second << std::endl; + LogTrace("VectorHitBuilderAlgorithm") << "A:\t upper local pos in the lower sof " << lparamsUpperInLower + << " with error: " << lparamsUpper.second << std::endl; + + bool ok = + checkClustersCompatibility(lparamsLower.first, lparamsUpper.first, lparamsLower.second, lparamsUpper.second); + + if (ok) { + AlgebraicSymMatrix22 covMat2Dzx; + double chi22Dzx = 0.0; + Local3DPoint pos2Dzx; + Local3DVector dir2Dzx; + fit2Dzx(lparamsLower.first, + lparamsUpperInLower, + lparamsLower.second, + lparamsUpper.second, + pos2Dzx, + dir2Dzx, + covMat2Dzx, + chi22Dzx); + LogTrace("VectorHitBuilderAlgorithm") << "\t pos2Dzx: " << pos2Dzx; + LogTrace("VectorHitBuilderAlgorithm") << "\t dir2Dzx: " << dir2Dzx; + LogTrace("VectorHitBuilderAlgorithm") << "\t cov2Dzx: " << covMat2Dzx; + VectorHit2D vh2Dzx = VectorHit2D(pos2Dzx, dir2Dzx, covMat2Dzx, chi22Dzx); + + AlgebraicSymMatrix22 covMat2Dzy; + double chi22Dzy = 0.0; + Local3DPoint pos2Dzy; + Local3DVector dir2Dzy; + fit2Dzy(lparamsLower.first, + lparamsUpperInLower, + lparamsLower.second, + lparamsUpper.second, + pos2Dzy, + dir2Dzy, + covMat2Dzy, + chi22Dzy); + LogTrace("VectorHitBuilderAlgorithm") << "\t pos2Dzy: " << pos2Dzy; + LogTrace("VectorHitBuilderAlgorithm") << "\t dir2Dzy: " << dir2Dzy; + LogTrace("VectorHitBuilderAlgorithm") << "\t cov2Dzy: " << covMat2Dzy; + VectorHit2D vh2Dzy = VectorHit2D(pos2Dzy, dir2Dzy, covMat2Dzy, chi22Dzy); + + OmniClusterRef lowerOmni(lower); + OmniClusterRef upperOmni(upper); + + Global3DPoint gPositionLower = VectorHit::phase2clusterGlobalPos(geomDetLower, lower); + Global3DPoint gPositionUpper = VectorHit::phase2clusterGlobalPos(geomDetUpper, upper); + GlobalError gErrorLower = VectorHit::phase2clusterGlobalPosErr(geomDetLower); + GlobalError gErrorUpper = VectorHit::phase2clusterGlobalPosErr(geomDetUpper); + + if (gPositionLower.perp() > gPositionUpper.perp()) { + std::swap(gPositionLower, gPositionUpper); + std::swap(gErrorLower, gErrorUpper); + } + + const CurvatureAndPhi curvatureAndPhi = curvatureANDphi(gPositionLower, gPositionUpper, gErrorLower, gErrorUpper); + VectorHit vh = VectorHit(*stack, + vh2Dzx, + vh2Dzy, + lowerOmni, + upperOmni, + curvatureAndPhi.curvature, + curvatureAndPhi.curvatureError, + curvatureAndPhi.phi); + return vh; + } + + return VectorHit(); +} + +void VectorHitBuilderAlgorithm::fit2Dzx(const Local3DPoint lpCI, + const Local3DPoint lpCO, + const LocalError leCI, + const LocalError leCO, + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const { + float x[2] = {lpCI.z(), lpCO.z()}; + float y[2] = {lpCI.x(), lpCO.x()}; + float sqCI = sqrt(leCI.xx()); + float sqCO = sqrt(leCO.xx()); + float sigy2[2] = {sqCI * sqCI, sqCO * sqCO}; + + fit(x, y, sigy2, pos, dir, covMatrix, chi2); + + return; +} + +void VectorHitBuilderAlgorithm::fit2Dzy(const Local3DPoint lpCI, + const Local3DPoint lpCO, + const LocalError leCI, + const LocalError leCO, + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const { + float x[2] = {lpCI.z(), lpCO.z()}; + float y[2] = {lpCI.y(), lpCO.y()}; + float sqCI = sqrt(leCI.yy()); + float sqCO = sqrt(leCO.yy()); + float sigy2[2] = {sqCI * sqCI, sqCO * sqCO}; + + fit(x, y, sigy2, pos, dir, covMatrix, chi2); + + return; +} + +void VectorHitBuilderAlgorithm::fit(float x[2], + float y[2], + float sigy2[2], + Local3DPoint& pos, + Local3DVector& dir, + AlgebraicSymMatrix22& covMatrix, + double& chi2) const { + float slope = 0.; + float intercept = 0.; + float covss = 0.; + float covii = 0.; + float covsi = 0.; + + linearFit(x, y, 2, sigy2, slope, intercept, covss, covii, covsi); + + covMatrix[0][0] = covss; // this is var(dy/dz) + covMatrix[1][1] = covii; // this is var(y) + covMatrix[1][0] = covsi; // this is cov(dy/dz,y) + + for (unsigned int j = 0; j < 2; j++) { + const double ypred = intercept + slope * x[j]; + const double dy = (y[j] - ypred) / sqrt(sigy2[j]); + chi2 += dy * dy; + } + + pos = Local3DPoint(intercept, 0., 0.); + //difference in z is the difference of the lowermost and the uppermost cluster z pos + float slopeZ = x[1] - x[0]; + dir = LocalVector(slope, 0., slopeZ); +} + +VectorHitBuilderAlgorithm::CurvatureAndPhi VectorHitBuilderAlgorithm::curvatureANDphi(Global3DPoint gPositionLower, + Global3DPoint gPositionUpper, + GlobalError gErrorLower, + GlobalError gErrorUpper) const { + VectorHitBuilderAlgorithm::CurvatureAndPhi result; + + float curvature = -999.; + float errorCurvature = -999.; + float phi = -999.; + + float h1 = gPositionLower.x() * gPositionUpper.y() - gPositionUpper.x() * gPositionLower.y(); + + //determine sign of curvature + AlgebraicVector2 n1; + n1[0] = -gPositionLower.y(); + n1[1] = gPositionLower.x(); + AlgebraicVector2 n2; + n2[0] = gPositionUpper.x() - gPositionLower.x(); + n2[1] = gPositionUpper.y() - gPositionLower.y(); + + double n3 = n1[0] * n2[0] + n1[1] * n2[1]; + double signCurv = -copysign(1.0, n3); + double phi1 = atan2(gPositionUpper.y() - gPositionLower.y(), gPositionUpper.x() - gPositionLower.x()); + + double x2Low = pow(gPositionLower.x(), 2); + double y2Low = pow(gPositionLower.y(), 2); + double x2Up = pow(gPositionUpper.x(), 2); + double y2Up = pow(gPositionUpper.y(), 2); + + if (h1 != 0) { + double h2 = 2 * h1; + double h2Inf = 1. / (2 * h1); + double r12 = gPositionLower.perp2(); + double r22 = gPositionUpper.perp2(); + double h3 = pow(n2[0], 2) + pow(n2[1], 2); + double h4 = -x2Low * gPositionUpper.x() + gPositionLower.x() * x2Up + gPositionLower.x() * y2Up - + gPositionUpper.x() * y2Low; + double h5 = + x2Low * gPositionUpper.y() - x2Up * gPositionLower.y() + y2Low * gPositionUpper.y() - gPositionLower.y() * y2Up; + + //radius of circle + double invRho2 = (4. * h1 * h1) / (r12 * r22 * h3); + curvature = sqrt(invRho2); + + //center of circle + double xcentre = h5 / h2; + double ycentre = h4 / h2; + + //to compute phi at the cluster points + double xtg = gPositionLower.y() - ycentre; + double ytg = -(gPositionLower.x() - xcentre); + + //to compute phi at the origin + phi = atan2(ytg, xtg); + + AlgebraicROOTObject<4, 4>::Matrix jacobian; + + double denom1 = 1. / sqrt(r12 * r22 * h3); + double denom2 = 1. / (pow(r12 * r22 * h3, 1.5)); + jacobian[0][0] = 1.0; // dx1/dx1 dx1/dy1 dx2/dx1 dy2/dx1 + jacobian[1][1] = 1.0; //dy1/dx1 dy1/dy1 dy2/dx1 dy2/dx1 + jacobian[2][0] = + -2. * ((h1 * (gPositionLower.x() * r22 * h3 + (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) * denom2 - + (gPositionUpper.y()) * denom1); // dkappa/dx1 + jacobian[2][1] = + -2. * ((gPositionUpper.x()) * denom1 + + (h1 * (gPositionLower.y() * r22 * h3 + r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) * + denom2); // dkappa/dy1 + jacobian[2][2] = + -2. * ((gPositionLower.y()) * denom1 + + (h1 * (gPositionUpper.x() * r12 * h3 - (gPositionLower.x() - gPositionUpper.x()) * r12 * r22)) * + denom2); // dkappa/dx2 + jacobian[2][3] = + -2. * ((h1 * (gPositionUpper.y() * r12 * h3 - r12 * r22 * (gPositionLower.y() - gPositionUpper.y()))) * denom2 - + (gPositionLower.x()) * denom1); // dkappa/dy2 + AlgebraicVector2 mVector; + //to compute phi at the cluster points + mVector[0] = (gPositionLower.y() - ycentre) * invRho2; // dphi/dxcentre + mVector[1] = -(gPositionLower.x() - xcentre) * invRho2; // dphi/dycentre + //to compute phi at the origin + + double h22Inv = 1. / pow(h2, 2); + + AlgebraicROOTObject<2, 4>::Matrix kMatrix; + kMatrix[0][0] = + 2. * ((gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionUpper.y() * h5) * h22Inv); // dxm/dx1 + kMatrix[0][1] = (2. * gPositionUpper.x() * h5) * h22Inv - + (x2Up + y2Up - 2. * gPositionLower.y() * gPositionUpper.y()) * h2Inf; // dxm/dy1 + kMatrix[0][2] = + 2. * ((gPositionLower.y() * h5) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf); // dxm/dx2 + kMatrix[0][3] = (x2Low + y2Low - 2. * gPositionUpper.y() * gPositionLower.y()) * h2Inf - + (2. * gPositionLower.x() * h5) * h22Inv; // dxm/dy2 + kMatrix[1][0] = (x2Up - 2. * gPositionLower.x() * gPositionUpper.x() + y2Up) * h2Inf - + (2. * gPositionUpper.y() * h4) * h22Inv; // dym/dx1 + kMatrix[1][1] = + 2. * ((gPositionUpper.x() * h4) * h22Inv - (gPositionUpper.x() * gPositionLower.y()) * h2Inf); // dym/dy1 + kMatrix[1][2] = (2. * gPositionLower.y() * h4) * h22Inv - + (x2Low - 2. * gPositionUpper.x() * gPositionLower.x() + y2Low) * h2Inf; // dym/dx2 + kMatrix[1][3] = + 2. * (gPositionLower.x() * gPositionUpper.y()) * h2Inf - (gPositionLower.x() * h4) * h22Inv; // dym/dy2 + + AlgebraicVector4 nMatrix = mVector * kMatrix; + jacobian[3][0] = nMatrix[0]; // dphi/(dx1,dy1,dx2,dy2) + jacobian[3][1] = nMatrix[1]; // dphi/(dx1,dy1,dx2,dy2) + jacobian[3][2] = nMatrix[2]; // dphi/(dx1,dy1,dx2,dy2) + jacobian[3][3] = nMatrix[3]; // dphi/(dx1,dy1,dx2,dy2) + + //assign correct sign to the curvature errors + if ((signCurv < 0 && curvature > 0) || (signCurv > 0 && curvature < 0)) { + curvature = -curvature; + for (int i = 0; i < 4; i++) { + jacobian[2][i] = -jacobian[2][i]; + } + } + + // bring phi in the same quadrant as phi1 + if (deltaPhi(phi, phi1) > M_PI / 2.) { + phi = phi + M_PI; + if (phi > M_PI) + phi = phi - 2. * M_PI; + } + + //computing the curvature error + AlgebraicVector4 curvatureJacobian; + for (int i = 0; i < 4; i++) { + curvatureJacobian[i] = jacobian[2][i]; + } + + AlgebraicROOTObject<4, 4>::Matrix gErrors; + + gErrors[0][0] = gErrorLower.cxx(); + gErrors[0][1] = gErrorLower.cyx(); + gErrors[1][0] = gErrorLower.cyx(); + gErrors[1][1] = gErrorLower.cyy(); + gErrors[2][2] = gErrorUpper.cxx(); + gErrors[2][3] = gErrorUpper.cyx(); + gErrors[3][2] = gErrorUpper.cyx(); + gErrors[3][3] = gErrorUpper.cyy(); + + AlgebraicVector4 temp = curvatureJacobian; + temp = temp * gErrors; + errorCurvature = temp[0] * curvatureJacobian[0] + temp[1] * curvatureJacobian[1] + temp[2] * curvatureJacobian[2] + + temp[3] * curvatureJacobian[3]; + } + + result.curvature = curvature; + result.curvatureError = errorCurvature; + result.phi = phi; + return result; +} diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithmBase.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithmBase.cc new file mode 100644 index 0000000000000..b875042c9fee1 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/src/VectorHitBuilderAlgorithmBase.cc @@ -0,0 +1,97 @@ +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithmBase.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "RecoLocalTracker/Records/interface/TkPhase2OTCPERecord.h" + +#include "Geometry/CommonTopologies/interface/PixelTopology.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +VectorHitBuilderAlgorithmBase::VectorHitBuilderAlgorithmBase( + const edm::ParameterSet& conf, + const TrackerGeometry* tkGeomProd, + const TrackerTopology* tkTopoProd, + const ClusterParameterEstimator* cpeProd) + : tkGeom_(tkGeomProd), + tkTopo_(tkTopoProd), + cpe_(cpeProd), + nMaxVHforeachStack_(conf.getParameter("maxVectorHitsInAStack")), + barrelCut_(conf.getParameter >("BarrelCut")), + endcapCut_(conf.getParameter >("EndcapCut")), + cpeTag_(conf.getParameter("CPE")) {} + +double VectorHitBuilderAlgorithmBase::computeParallaxCorrection(const PixelGeomDetUnit* geomDetUnit_low, + const Point3DBase& lPosClu_low, + const PixelGeomDetUnit* geomDetUnit_upp, + const Point3DBase& lPosClu_upp) const { + double parallCorr = 0.0; + Global3DPoint origin(0, 0, 0); + Global3DPoint gPosClu_low = geomDetUnit_low->surface().toGlobal(lPosClu_low); + GlobalVector gV = gPosClu_low - origin; + LogTrace("VectorHitsBuilderValidation") << " global vector passing to the origin:" << gV; + + LocalVector lV = geomDetUnit_low->surface().toLocal(gV); + LogTrace("VectorHitsBuilderValidation") + << " local vector passing to the origin (in the lower detector system of reference):" << lV; + LocalVector lV_norm = lV / lV.z(); + LogTrace("VectorHitsBuilderValidation") + << " normalized local vector passing to the origin (in low the lower detector system of reference):" << lV_norm; + + Global3DPoint gPosClu_upp = geomDetUnit_upp->surface().toGlobal(lPosClu_upp); + Local3DPoint lPosClu_uppInLow = geomDetUnit_low->surface().toLocal(gPosClu_upp); + parallCorr = lV_norm.x() * lPosClu_uppInLow.z(); + + return parallCorr; +} + +void VectorHitBuilderAlgorithmBase::printClusters(const edmNew::DetSetVector& clusters) const { + int nCluster = 0; + int numberOfDSV = 0; + for (const auto& DSViter : clusters) { + ++numberOfDSV; + // Loop over the clusters in the detector unit + for (const auto& clustIt : DSViter) { + nCluster++; + // get the detector unit's id + const GeomDetUnit* geomDetUnit(tkGeom_->idToDetUnit(DSViter.detId())); + if (!geomDetUnit) + return; + printCluster(geomDetUnit, &clustIt); + } + } + LogDebug("VectorHitBuilder") << " Number of input clusters: " << nCluster << std::endl; +} + +void VectorHitBuilderAlgorithmBase::printCluster(const GeomDet* geomDetUnit, + const Phase2TrackerCluster1D* clustIt) const { + if (!geomDetUnit) + return; + const PixelGeomDetUnit* pixelGeomDetUnit = dynamic_cast(geomDetUnit); + const PixelTopology& topol = pixelGeomDetUnit->specificTopology(); + if (!pixelGeomDetUnit) + return; + + unsigned int layer = tkTopo_->layer(geomDetUnit->geographicalId()); + unsigned int module = tkTopo_->module(geomDetUnit->geographicalId()); + LogTrace("VectorHitBuilder") << "Layer:" << layer << " and DetId: " << geomDetUnit->geographicalId().rawId() + << std::endl; + TrackerGeometry::ModuleType mType = tkGeom_->getDetectorType(geomDetUnit->geographicalId()); + if (mType == TrackerGeometry::ModuleType::Ph2PSP) + LogTrace("VectorHitBuilder") << "Pixel cluster (module:" << module << ") " << std::endl; + else if (mType == TrackerGeometry::ModuleType::Ph2SS || mType == TrackerGeometry::ModuleType::Ph2PSS) + LogTrace("VectorHitBuilder") << "Strip cluster (module:" << module << ") " << std::endl; + else + LogTrace("VectorHitBuilder") << "no module?!" << std::endl; + LogTrace("VectorHitBuilder") << "with pitch:" << topol.pitch().first << " , " << topol.pitch().second << std::endl; + LogTrace("VectorHitBuilder") << " and width:" << pixelGeomDetUnit->surface().bounds().width() + << " , lenght:" << pixelGeomDetUnit->surface().bounds().length() << std::endl; + + auto&& lparams = cpe_->localParameters(*clustIt, *pixelGeomDetUnit); + Global3DPoint gparams = pixelGeomDetUnit->surface().toGlobal(lparams.first); + + LogTrace("VectorHitBuilder") << "\t global pos " << gparams << std::endl; + LogTrace("VectorHitBuilder") << "\t local pos " << lparams.first << "with err " << lparams.second << std::endl; + LogTrace("VectorHitBuilder") << std::endl; + + return; +} diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/BuildFile.xml b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/BuildFile.xml new file mode 100644 index 0000000000000..847e241c1ee43 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/BuildFile.xml @@ -0,0 +1,30 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/ClustersValidationTGraph.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/ClustersValidationTGraph.cc new file mode 100644 index 0000000000000..bb12a15caee68 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/ClustersValidationTGraph.cc @@ -0,0 +1,551 @@ +#include +#include +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h" +#include "DataFormats/SiPixelDigi/interface/PixelDigi.h" + +#include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" + +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "CommonTools/Utils/interface/TFileDirectory.h" + +#include +#include +#include +#include + +struct ClusterHistos { + THStack* numberClustersMixed; + TH1F* numberClusterPixel; + TH1F* numberClusterStrip; + + THStack* clustersSizeMixed; + TH1F* clusterSizePixel; + TH1F* clusterSizeStrip; + + TGraph* globalPosXY[3]; + TGraph* localPosXY[3]; + + TH1F* deltaXClusterSimHits[3]; + TH1F* deltaYClusterSimHits[3]; + + TH1F* deltaXClusterSimHits_P[3]; + TH1F* deltaYClusterSimHits_P[3]; + + TH1F* digiEfficiency[3]; + + TH1F* primarySimHits; + TH1F* otherSimHits; +}; + +class Phase2TrackerClusterizerValidationTGraph : public edm::one::EDAnalyzer { +public: + typedef std::map > SimHitsMap; + typedef std::map SimTracksMap; + + explicit Phase2TrackerClusterizerValidationTGraph(const edm::ParameterSet&); + ~Phase2TrackerClusterizerValidationTGraph(); + void beginJob(); + void endJob(); + void analyze(const edm::Event&, const edm::EventSetup&); + static void fillDescriptions(edm::ConfigurationDescriptions&); + +private: + std::map::iterator createLayerHistograms(unsigned int); + unsigned int getLayerNumber(const DetId&, const TrackerTopology*); + unsigned int getModuleNumber(const DetId&, const TrackerTopology*); + unsigned int getSimTrackId(const edm::Handle >&, const DetId&, unsigned int); + edm::EDGetTokenT > srcClu_; + edm::EDGetTokenT > siphase2OTSimLinksToken_; + edm::EDGetTokenT simHitsToken_; + edm::EDGetTokenT simTracksToken_; + edm::EDGetTokenT simVerticesToken_; + const TrackerGeometry* tkGeom_; + const TrackerTopology* tkTopo_; + + TTree* tree_; + TGraph* trackerLayout_[3]; + TGraph* trackerLayoutXY_[3]; + TGraph* trackerLayoutXYBar_; + TGraph* trackerLayoutXYEC_; + + std::map histograms_; +}; + +Phase2TrackerClusterizerValidationTGraph::Phase2TrackerClusterizerValidationTGraph(const edm::ParameterSet& conf) { + srcClu_ = + consumes >(edm::InputTag(conf.getParameter("src"))); + siphase2OTSimLinksToken_ = consumes >(conf.getParameter("links")); + simHitsToken_ = consumes(edm::InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof")); + simTracksToken_ = consumes(edm::InputTag("g4SimHits")); + simVerticesToken_ = consumes(edm::InputTag("g4SimHits")); +} + +Phase2TrackerClusterizerValidationTGraph::~Phase2TrackerClusterizerValidationTGraph() {} + +void Phase2TrackerClusterizerValidationTGraph::beginJob() { + edm::Service fs; + fs->file().cd("/"); + TFileDirectory td = fs->mkdir("Common"); + //Create common ntuple + tree_ = td.make("Phase2TrackerClusters", "Phase2TrackerClusters"); + // Create common histograms + trackerLayout_[0] = td.make(); + trackerLayout_[0]->SetName("RVsZ_Mixed"); + trackerLayout_[1] = td.make(); + trackerLayout_[1]->SetName("RVsZ_Pixel"); + trackerLayout_[2] = td.make(); + trackerLayout_[2]->SetName("RVsZ_Strip"); + trackerLayoutXY_[0] = td.make(); + trackerLayoutXY_[0]->SetName("YVsX_Mixed"); + trackerLayoutXY_[1] = td.make(); + trackerLayoutXY_[1]->SetName("YVsX_Pixel"); + trackerLayoutXY_[2] = td.make(); + trackerLayoutXY_[2]->SetName("YVsX_Strip"); + trackerLayoutXYBar_ = td.make(); + trackerLayoutXYBar_->SetName("YVsXBar"); + trackerLayoutXYEC_ = td.make(); + trackerLayoutXYEC_->SetName("YVsXEC"); +} + +void Phase2TrackerClusterizerValidationTGraph::endJob() {} + +void Phase2TrackerClusterizerValidationTGraph::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) { + // Get the needed objects + + // Get the clusters + edm::Handle clusters; + event.getByToken(srcClu_, clusters); + + // Get the Phase2 DigiSimLink + edm::Handle > siphase2SimLinks; + event.getByToken(siphase2OTSimLinksToken_, siphase2SimLinks); + + // Get the SimHits + edm::Handle simHitsRaw; + event.getByToken(simHitsToken_, simHitsRaw); + + // Get the SimTracks + edm::Handle simTracksRaw; + event.getByToken(simTracksToken_, simTracksRaw); + + // Get the SimVertex + edm::Handle simVertices; + event.getByToken(simVerticesToken_, simVertices); + + // Get the geometry + edm::ESHandle geomHandle; + eventSetup.get().get(geomHandle); + tkGeom_ = &(*geomHandle); + edm::ESHandle tTopoHandle; + eventSetup.get().get(tTopoHandle); + tkTopo_ = tTopoHandle.product(); + + //set up for tree + int layer_number; + int module_id; + int module_number; + int module_type; //1: pixel, 2: strip + float x_global, y_global, z_global; + float x_local, y_local, z_local; + + tree_->Branch("layer_number", &layer_number, "layer_number/I"); + tree_->Branch("module_id", &module_id, "module_id/I"); + tree_->Branch("module_type", &module_type, "module_type/I"); + tree_->Branch("module_number", &module_number, "module_number/I"); + tree_->Branch("x_global", &x_global, "x_global/F"); + tree_->Branch("y_global", &y_global, "y_global/F"); + tree_->Branch("z_global", &z_global, "z_global/F"); + tree_->Branch("x_local", &x_local, "x_local/F"); + tree_->Branch("y_local", &y_local, "y_local/F"); + tree_->Branch("z_local", &z_local, "z_local/F"); + + // Rearrange the simTracks for ease of use + SimTracksMap simTracks; + for (const auto& simTrackIt : *simTracksRaw) + simTracks.emplace(std::pair(simTrackIt.trackId(), simTrackIt)); + + // Rearrange the simHits by detUnit + + // Rearrange the simHits for ease of use + SimHitsMap simHitsDetUnit; + SimHitsMap simHitsTrackId; + for (const auto& simHitIt : *simHitsRaw) { + SimHitsMap::iterator simHitsDetUnitIt(simHitsDetUnit.find(simHitIt.detUnitId())); + if (simHitsDetUnitIt == simHitsDetUnit.end()) { + std::pair newIt(simHitsDetUnit.insert( + std::pair >(simHitIt.detUnitId(), std::vector()))); + simHitsDetUnitIt = newIt.first; + } + simHitsDetUnitIt->second.push_back(simHitIt); + SimHitsMap::iterator simHitsTrackIdIt(simHitsTrackId.find(simHitIt.trackId())); + if (simHitsTrackIdIt == simHitsTrackId.end()) { + std::pair newIt(simHitsTrackId.insert( + std::pair >(simHitIt.trackId(), std::vector()))); + simHitsTrackIdIt = newIt.first; + } + simHitsTrackIdIt->second.push_back(simHitIt); + } + + // ValidationTGraph + unsigned int nClustersTot(0), nClustersPixelTot(0), nClustersStripTot(0); + + // Loop over modules + for (const auto& DSViter : *clusters) { + // Get the detector unit's id + unsigned int rawid(DSViter.detId()); + module_id = rawid; + DetId detId(rawid); + + layer_number = getLayerNumber(detId, tkTopo_); + module_number = getModuleNumber(detId, tkTopo_); + unsigned int layer(getLayerNumber(detId, tkTopo_)); + + // Get the geometry of the tracker + const GeomDetUnit* geomDetUnit(tkGeom_->idToDetUnit(detId)); + const PixelGeomDetUnit* theGeomDet = dynamic_cast(geomDetUnit); + const PixelTopology& topol = theGeomDet->specificTopology(); + + if (!geomDetUnit) + break; + + // Create histograms for the layer if they do not yet exist + std::map::iterator histogramLayer(histograms_.find(layer)); + if (histogramLayer == histograms_.end()) + histogramLayer = createLayerHistograms(layer); + + // Number of clusters + unsigned int nClustersPixel(0), nClustersStrip(0); + + // Loop over the clusters in the detector unit + for (const auto& clustIt : DSViter) { + // Cluster related variables + MeasurementPoint mpClu(clustIt.center(), clustIt.column() + 0.5); + Local3DPoint localPosClu = geomDetUnit->topology().localPosition(mpClu); + x_local = localPosClu.x(); + y_local = localPosClu.y(); + z_local = localPosClu.z(); + + Global3DPoint globalPosClu = geomDetUnit->surface().toGlobal(localPosClu); + x_global = globalPosClu.x(); + y_global = globalPosClu.y(); + z_global = globalPosClu.z(); + + // Fill the position histograms + trackerLayout_[0]->SetPoint(nClustersTot, globalPosClu.z(), globalPosClu.perp()); + trackerLayoutXY_[0]->SetPoint(nClustersTot, globalPosClu.x(), globalPosClu.y()); + + if (layer < 100) + trackerLayoutXYBar_->SetPoint(nClustersTot, globalPosClu.x(), globalPosClu.y()); + else + trackerLayoutXYEC_->SetPoint(nClustersTot, globalPosClu.x(), globalPosClu.y()); + + histogramLayer->second.localPosXY[0]->SetPoint(nClustersTot, localPosClu.x(), localPosClu.y()); + histogramLayer->second.globalPosXY[0]->SetPoint(nClustersTot, globalPosClu.x(), globalPosClu.y()); + + // Pixel module + if (topol.ncolumns() == 32) { + module_type = 1; + trackerLayout_[1]->SetPoint(nClustersPixelTot, globalPosClu.z(), globalPosClu.perp()); + trackerLayoutXY_[1]->SetPoint(nClustersPixelTot, globalPosClu.x(), globalPosClu.y()); + + histogramLayer->second.localPosXY[1]->SetPoint(nClustersPixelTot, localPosClu.x(), localPosClu.y()); + histogramLayer->second.globalPosXY[1]->SetPoint(nClustersPixelTot, globalPosClu.x(), globalPosClu.y()); + histogramLayer->second.clusterSizePixel->Fill(clustIt.size()); + ++nClustersPixel; + ++nClustersPixelTot; + } + // Strip module + else if (topol.ncolumns() == 2) { + module_type = 2; + trackerLayout_[2]->SetPoint(nClustersStripTot, globalPosClu.z(), globalPosClu.perp()); + trackerLayoutXY_[2]->SetPoint(nClustersStripTot, globalPosClu.x(), globalPosClu.y()); + + histogramLayer->second.localPosXY[2]->SetPoint(nClustersStripTot, localPosClu.x(), localPosClu.y()); + histogramLayer->second.globalPosXY[2]->SetPoint(nClustersStripTot, globalPosClu.x(), globalPosClu.y()); + histogramLayer->second.clusterSizeStrip->Fill(clustIt.size()); + ++nClustersStrip; + ++nClustersStripTot; + } + + // * Digis related variables + + std::vector clusterSimTrackIds; + + // Get all the simTracks that form the cluster + for (unsigned int i(0); i < clustIt.size(); ++i) { + unsigned int channel(PixelDigi::pixelToChannel( + clustIt.firstRow() + i, + clustIt + .column())); // Here we have to use the old pixelToChannel function (not Phase2TrackerDigi but PixelDigi), change this when using new Digis + unsigned int simTrackId(getSimTrackId(siphase2SimLinks, detId, channel)); + clusterSimTrackIds.push_back(simTrackId); + } + } + + if (nClustersPixel) + histogramLayer->second.numberClusterPixel->Fill(nClustersPixel); + if (nClustersStrip) + histogramLayer->second.numberClusterStrip->Fill(nClustersStrip); + nClustersTot++; + tree_->Fill(); + } +} + +// Create the histograms +std::map::iterator Phase2TrackerClusterizerValidationTGraph::createLayerHistograms( + unsigned int ival) { + std::ostringstream fname1, fname2; + + edm::Service fs; + fs->file().cd("/"); + + std::string tag; + unsigned int id; + if (ival < 100) { + id = ival; + fname1 << "Barrel"; + fname2 << "Layer_" << id; + tag = "_layer_"; + } else { + int side = ival / 100; + id = ival - side * 100; + fname1 << "EndCap_Side_" << side; + fname2 << "Disc_" << id; + tag = "_disc_"; + } + + TFileDirectory td1 = fs->mkdir(fname1.str().c_str()); + TFileDirectory td = td1.mkdir(fname2.str().c_str()); + + ClusterHistos local_histos; + + std::ostringstream histoName; + + /* + * Number of clusters + */ + + histoName.str(""); + histoName << "Number_Clusters_Pixel" << tag.c_str() << id; + local_histos.numberClusterPixel = td.make(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.); + local_histos.numberClusterPixel->SetFillColor(kAzure + 7); + + histoName.str(""); + histoName << "Number_Clusters_Strip" << tag.c_str() << id; + local_histos.numberClusterStrip = td.make(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.); + local_histos.numberClusterStrip->SetFillColor(kOrange - 3); + + histoName.str(""); + histoName << "Number_Clusters_Mixed" << tag.c_str() << id; + local_histos.numberClustersMixed = td.make(histoName.str().c_str(), histoName.str().c_str()); + local_histos.numberClustersMixed->Add(local_histos.numberClusterPixel); + local_histos.numberClustersMixed->Add(local_histos.numberClusterStrip); + + /* + * Cluster size + */ + + histoName.str(""); + histoName << "Cluster_Size_Pixel" << tag.c_str() << id; + local_histos.clusterSizePixel = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + local_histos.clusterSizePixel->SetFillColor(kAzure + 7); + + histoName.str(""); + histoName << "Cluster_Size_Strip" << tag.c_str() << id; + local_histos.clusterSizeStrip = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + local_histos.clusterSizeStrip->SetFillColor(kOrange - 3); + + histoName.str(""); + histoName << "Cluster_Size_Mixed" << tag.c_str() << id; + local_histos.clustersSizeMixed = td.make(histoName.str().c_str(), histoName.str().c_str()); + local_histos.clustersSizeMixed->Add(local_histos.clusterSizePixel); + local_histos.clustersSizeMixed->Add(local_histos.clusterSizeStrip); + + /* + * Local and Global positions + */ + + histoName.str(""); + histoName << "Local_Position_XY_Mixed" << tag.c_str() << id; + local_histos.localPosXY[0] = td.make(); + local_histos.localPosXY[0]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Local_Position_XY_Pixel" << tag.c_str() << id; + local_histos.localPosXY[1] = td.make(); + local_histos.localPosXY[1]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Local_Position_XY_Strip" << tag.c_str() << id; + local_histos.localPosXY[2] = td.make(); + local_histos.localPosXY[2]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_Mixed" << tag.c_str() << id; + local_histos.globalPosXY[0] = td.make(); + local_histos.globalPosXY[0]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_Pixel" << tag.c_str() << id; + local_histos.globalPosXY[1] = td.make(); + local_histos.globalPosXY[1]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_Strip" << tag.c_str() << id; + local_histos.globalPosXY[2] = td.make(); + local_histos.globalPosXY[2]->SetName(histoName.str().c_str()); + + /* + * Delta positions with SimHits + */ + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Mixed" << tag.c_str() << id; + local_histos.deltaXClusterSimHits[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Pixel" << tag.c_str() << id; + local_histos.deltaXClusterSimHits[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Strip" << tag.c_str() << id; + local_histos.deltaXClusterSimHits[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Mixed" << tag.c_str() << id; + local_histos.deltaYClusterSimHits[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Pixel" << tag.c_str() << id; + local_histos.deltaYClusterSimHits[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Strip" << tag.c_str() << id; + local_histos.deltaYClusterSimHits[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + /* + * Delta position with simHits for primary tracks only + */ + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Mixed_P" << tag.c_str() << id; + local_histos.deltaXClusterSimHits_P[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Pixel_P" << tag.c_str() << id; + local_histos.deltaXClusterSimHits_P[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_Cluster_SimHits_Strip_P" << tag.c_str() << id; + local_histos.deltaXClusterSimHits_P[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Mixed_P" << tag.c_str() << id; + local_histos.deltaYClusterSimHits_P[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Pixel_P" << tag.c_str() << id; + local_histos.deltaYClusterSimHits_P[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_Cluster_SimHits_Strip_P" << tag.c_str() << id; + local_histos.deltaYClusterSimHits_P[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + /* + * Information on the Digis per cluster + */ + + histoName.str(""); + histoName << "Primary_Digis" << tag.c_str() << id; + local_histos.primarySimHits = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + + histoName.str(""); + histoName << "Other_Digis" << tag.c_str() << id; + local_histos.otherSimHits = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + + /* + * End + */ + + std::pair::iterator, bool> insertedIt( + histograms_.insert(std::make_pair(ival, local_histos))); + fs->file().cd("/"); + + return insertedIt.first; +} + +unsigned int Phase2TrackerClusterizerValidationTGraph::getLayerNumber(const DetId& detid, const TrackerTopology* topo) { + if (detid.det() == DetId::Tracker) { + if (detid.subdetId() == PixelSubdetector::PixelBarrel) + return (topo->pxbLayer(detid)); + else if (detid.subdetId() == PixelSubdetector::PixelEndcap) + return (100 * topo->pxfSide(detid) + topo->pxfDisk(detid)); + } + return 999; +} + +unsigned int Phase2TrackerClusterizerValidationTGraph::getModuleNumber(const DetId& detid, + const TrackerTopology* topo) { + if (detid.det() == DetId::Tracker) { + if (detid.subdetId() == PixelSubdetector::PixelBarrel) { + return (topo->pxbModule(detid)); + } else if (detid.subdetId() == PixelSubdetector::PixelEndcap) { + return (topo->pxfModule(detid)); + } else + return 999; + } + return 999; +} + +unsigned int Phase2TrackerClusterizerValidationTGraph::getSimTrackId( + const edm::Handle >& siphase2SimLinks, + const DetId& detId, + unsigned int channel) { + edm::DetSetVector::const_iterator DSViter(siphase2SimLinks->find(detId)); + if (DSViter == siphase2SimLinks->end()) + return 0; + for (const auto& it : DSViter->data) { + if (channel == it.channel()) + return it.SimTrackId(); + } + return 0; +} + +void Phase2TrackerClusterizerValidationTGraph::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("src", "siPhase2Clusters"); + desc.add("links", edm::InputTag("simSiPixelDigis", "Tracker")); + descriptions.add("phase2TrackerClusterizerValidationTGraph", desc); +} + +DEFINE_FWK_MODULE(Phase2TrackerClusterizerValidationTGraph); diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/Clusters_productionAndTesting.py b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/Clusters_productionAndTesting.py new file mode 100644 index 0000000000000..cac36312ac6a7 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/Clusters_productionAndTesting.py @@ -0,0 +1,119 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('RECO',eras.Phase2C9) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +#process.load('Configuration.StandardSequences.Validation_cff') +#process.load('DQMOffline.Configuration.DQMOfflineMC_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +#adding only recolocalreco +process.load('RecoLocalTracker.Configuration.RecoLocalTracker_cff') + +# import VectorHitBuilder +process.load('RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2VectorHits_cfi') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/863B48BB-03ED-0548-AA60-1269291ED1E6.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/9B6E3A66-B330-8E42-B85E-96A9952A002E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/29706FE9-16C9-CE4F-B744-66E07B250D1E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AAB16BEE-B0CE-644A-8E96-35236D793C04.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/36F89E35-DE2F-174B-95B7-7A9423DED2D8.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/1BE0A565-4F64-8D42-A5D5-62692F64F0A5.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/8671B5C9-DE1F-8B4C-8A1C-9C99898FE191.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/6CA0E490-73F4-1147-906D-050B8B3A3134.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AFAB961F-E31F-064D-9031-BEAA10702345.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/2AB6718E-EEA6-494B-AC25-B59CF36DF941.root', + ), + secondaryFileNames = cms.untracked.vstring(), + skipEvents = cms.untracked.uint32(0) +) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('step3 nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOSIMoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('GEN-SIM-RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('file:step3_1event.root'), + outputCommands = process.RECOSIMEventContent.outputCommands, + #outputCommands = cms.untracked.vstring( ('keep *') ), + splitLevel = cms.untracked.int32(0) +) + +# debug +process.MessageLogger = cms.Service('MessageLogger', + debugModules = cms.untracked.vstring('siPhase2Clusters'), + destinations = cms.untracked.vstring('cout'), + cout = cms.untracked.PSet( + threshold = cms.untracked.string('ERROR') + ) +) + +# Analyzer +# Analyzer +process.analysis = cms.EDAnalyzer('Phase2TrackerClusterizerValidation', + src = cms.InputTag("siPhase2Clusters"), + links = cms.InputTag("simSiPixelDigis", "Tracker"), + simhitsbarrel = cms.InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof"), + simhitsendcap = cms.InputTag("g4SimHits", "TrackerHitsPixelEndcapLowTof"), + simtracks = cms.InputTag("g4SimHits"), + ECasRings = cms.bool(True), + SimTrackMinPt = cms.double(2.) +) + +#process.analysis = cms.EDAnalyzer('Phase2TrackerClusterizerValidationTGraph', +# src = cms.string("siPhase2Clusters"), +# links = cms.InputTag("simSiPixelDigis", "Tracker") +#) +process.TFileService = cms.Service('TFileService', + fileName = cms.string('file:Clusters_validation.root') +) + + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.trackerlocalreco_step = cms.Path(process.trackerlocalreco) +process.analysis_step = cms.Path(process.analysis) +process.RECOSIMoutput_step = cms.EndPath(process.RECOSIMoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.trackerlocalreco_step,process.RECOSIMoutput_step, process.analysis_step) + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_SeedingOT_productionAndTesting.py b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_SeedingOT_productionAndTesting.py new file mode 100644 index 0000000000000..ac5fdd513eb95 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_SeedingOT_productionAndTesting.py @@ -0,0 +1,121 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('RECO',eras.Phase2C9) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +#adding only recolocalreco +process.load('RecoLocalTracker.Configuration.RecoLocalTracker_cff') + +# import VectorHitBuilder +process.load('RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2VectorHits_cfi') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/863B48BB-03ED-0548-AA60-1269291ED1E6.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/9B6E3A66-B330-8E42-B85E-96A9952A002E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/29706FE9-16C9-CE4F-B744-66E07B250D1E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AAB16BEE-B0CE-644A-8E96-35236D793C04.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/36F89E35-DE2F-174B-95B7-7A9423DED2D8.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/1BE0A565-4F64-8D42-A5D5-62692F64F0A5.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/8671B5C9-DE1F-8B4C-8A1C-9C99898FE191.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/6CA0E490-73F4-1147-906D-050B8B3A3134.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AFAB961F-E31F-064D-9031-BEAA10702345.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/2AB6718E-EEA6-494B-AC25-B59CF36DF941.root', + ), + secondaryFileNames = cms.untracked.vstring(), + skipEvents = cms.untracked.uint32(0) +) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('step3 nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOSIMoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('GEN-SIM-RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('file:step3_1event.root'), + outputCommands = process.RECOSIMEventContent.outputCommands, + splitLevel = cms.untracked.int32(0) +) + +# Analyzer +process.load('RecoLocalTracker.SiPhase2VectorHitBuilder.vectorHitsBuilderValidation_cfi') + +process.TFileService = cms.Service('TFileService', + fileName = cms.string('file:vh_validation_tilted.root') +) + +from RecoTracker.TkSeedGenerator.SeedingOTEDProducer_cfi import SeedingOTEDProducer as _SeedingOTEDProducer +process.phase2SeedingOTEDProducer = _SeedingOTEDProducer.clone() +process.initialStepSeeds = _SeedingOTEDProducer.clone() + +process.load('RecoLocalTracker.Phase2TrackerRecHits.Phase2StripCPEGeometricESProducer_cfi') + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + + +# debug +#process.MessageLogger = cms.Service("MessageLogger", +# destinations = cms.untracked.vstring("debugVH_tilted"), +# debugModules = cms.untracked.vstring("*"), +# categories = cms.untracked.vstring("VectorHitBuilderEDProducer","VectorHitBuilderAlgorithm","VectorHitsBuilderValidation","CkfPattern"), +# debugVH_tilted = cms.untracked.PSet(threshold = cms.untracked.string("DEBUG"), +# DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(0)), +# default = cms.untracked.PSet(limit = cms.untracked.int32(0)), +# VectorHitBuilderEDProducer = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# VectorHitBuilderAlgorithm = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# CkfPattern = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# VectorHitsBuilderValidation = cms.untracked.PSet(limit = cms.untracked.int32(-1)) +# ) +# ) + + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.trackerlocalreco_step = cms.Path(process.trackerlocalreco+process.siPixelClusters+process.siPhase2VectorHits) +process.seedingOT_step = cms.Path(process.MeasurementTrackerEvent+process.offlineBeamSpot+process.phase2SeedingOTEDProducer) +process.analysis_step = cms.Path(process.vectorHitsBuilderValidation) +process.RECOSIMoutput_step = cms.EndPath(process.RECOSIMoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.trackerlocalreco_step,process.seedingOT_step,process.RECOSIMoutput_step, process.analysis_step) + +# customisation of the process. + + +# End of customisation functions + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_combinatorialStudies_PU200.py b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_combinatorialStudies_PU200.py new file mode 100644 index 0000000000000..4a6c601fa2ecd --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_combinatorialStudies_PU200.py @@ -0,0 +1,135 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('RECO',eras.Phase2C9) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +#process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('SimGeneral.MixingModule.mix_POISSON_average_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +#adding only recolocalreco +process.load('RecoLocalTracker.Configuration.RecoLocalTracker_cff') + +# import VectorHitBuilder +process.load('RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2VectorHits_cfi') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/C7A95D39-8FAB-D746-8C06-B11AECA047DD.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/E4D19FDC-9192-2F44-A4FA-21AB5A93344C.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/04E557E8-C456-B148-9EC6-2E5E40BD03C6.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/284ACE29-A506-CD40-9059-899981D7931C.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/5ACFC231-EBF6-2B4F-A88F-D01171F506CE.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/8E76FEBA-FF9B-544D-9C17-C5B31777FE39.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/2958EE82-39A5-154A-9650-5E97AAD44B3E.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/943CFE41-83ED-F04B-AA74-435821572336.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/01D965C0-05A4-C543-A676-7C6C67DC2175.root', +'/store/relval/CMSSW_11_2_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/20000/7D171FAD-4E5B-824A-A0DF-07175F0A9550.root', +), + #fileNames = cms.untracked.vstring('file:step2.root'), + secondaryFileNames = cms.untracked.vstring(), + skipEvents = cms.untracked.uint32(0) +) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('step3 nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOSIMoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('GEN-SIM-RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('file:step3_PU200.root'), + outputCommands = process.RECOSIMEventContent.outputCommands, + #outputCommands = cms.untracked.vstring( ('keep *') ), + splitLevel = cms.untracked.int32(0) +) + +# debug +process.MessageLogger = cms.Service("MessageLogger", + destinations = cms.untracked.vstring("debugVH_PU200"), + debugModules = cms.untracked.vstring("*"), + categories = cms.untracked.vstring("VectorHitsBuilderValidation"), + debugVH_PU200 = cms.untracked.PSet(threshold = cms.untracked.string("DEBUG"), + DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(0)), + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), + VectorHitsBuilderValidation = cms.untracked.PSet(limit = cms.untracked.int32(-1)) + ) + ) + +# Analyzer +process.analysis = cms.EDAnalyzer('VectorHitsBuilderValidation', + src = cms.string("siPhase2Clusters"), + VH_acc = cms.InputTag("siPhase2VectorHits", "accepted"), + VH_rej = cms.InputTag("siPhase2VectorHits", "rejected"), + CPE = cms.ESInputTag("phase2StripCPEESProducer", "Phase2StripCPE"), + links = cms.InputTag("simSiPixelDigis", "Tracker"), + trackingParticleSrc = cms.InputTag('mix', 'MergedTrackTruth'), +) +process.TFileService = cms.Service('TFileService', + fileName = cms.string('file:VHs_validation_PU200_new.root') +) + + +# Other statements +process.mix.input.nbPileupEvents.averageNumber = cms.double(200.000000) +process.mix.bunchspace = cms.int32(25) +process.mix.minBunch = cms.int32(-3) +process.mix.maxBunch = cms.int32(3) +process.mix.input.fileNames = cms.untracked.vstring([ +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/09AA31DF-737A-B142-A82B-6E48B4B965CD.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/A18897B2-7509-C845-B83A-CEE4AC55FF86.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/184F91C7-EA57-2B49-9E85-D1F494152CDD.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/A505AC6C-6DEA-484C-8FFC-1EA4B54892F5.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/5636AD91-D035-0B44-8D53-2B4FFE6A7722.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/03D1F17E-33EE-A94C-A869-00071F1F45C7.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/EB99E4B3-838A-CC4C-B1B2-8A9D1E1EB453.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/7F44A9B9-38C1-484F-9B4B-DCD0FD7224A8.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/2A3CCF08-CDD1-6246-827F-27CD0CC42478.root', +'/store/relval/CMSSW_11_2_0_pre5/RelValMinBias_14TeV/GEN-SIM-DIGI-RAW/110X_mcRun4_realistic_v3_2026D49noPU-v1/20000/0745B823-5DB0-B74F-B0A6-78E24ED39EE5.root', +]) +process.mix.playback = True +process.mix.digitizers = cms.PSet() +for a in process.aliases: delattr(process, a) +process.RandomNumberGeneratorService.restoreStateLabel=cms.untracked.string("randomEngineStateProducer") +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.trackerlocalreco_step = cms.Path(process.trackerlocalreco+process.siPhase2VectorHits) +process.analysis_step = cms.Path(process.analysis) +process.RECOSIMoutput_step = cms.EndPath(process.RECOSIMoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.trackerlocalreco_step,process.RECOSIMoutput_step, process.analysis_step) + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_productionAndTesting.py b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_productionAndTesting.py new file mode 100644 index 0000000000000..cbb75e27d744a --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VHs_productionAndTesting.py @@ -0,0 +1,114 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('RECO',eras.Phase2C9) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +#adding only recolocalreco +process.load('RecoLocalTracker.Configuration.RecoLocalTracker_cff') + +# import VectorHitBuilder +process.load('RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2VectorHits_cfi') + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/863B48BB-03ED-0548-AA60-1269291ED1E6.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/9B6E3A66-B330-8E42-B85E-96A9952A002E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/29706FE9-16C9-CE4F-B744-66E07B250D1E.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AAB16BEE-B0CE-644A-8E96-35236D793C04.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/36F89E35-DE2F-174B-95B7-7A9423DED2D8.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/1BE0A565-4F64-8D42-A5D5-62692F64F0A5.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/8671B5C9-DE1F-8B4C-8A1C-9C99898FE191.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/6CA0E490-73F4-1147-906D-050B8B3A3134.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/AFAB961F-E31F-064D-9031-BEAA10702345.root', + '/store/relval/CMSSW_11_2_0_pre3/RelValSingleMuFlatPt2To100/GEN-SIM-DIGI-RAW/PU25ns_110X_mcRun4_realistic_v3_2026D49PU200-v1/20000/2AB6718E-EEA6-494B-AC25-B59CF36DF941.root', + ), + secondaryFileNames = cms.untracked.vstring(), + skipEvents = cms.untracked.uint32(0) +) + +process.options = cms.untracked.PSet( + +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('step3 nevts:10'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOSIMoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('GEN-SIM-RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('file:step3_1event.root'), + outputCommands = process.RECOSIMEventContent.outputCommands, + splitLevel = cms.untracked.int32(0) +) + +# debug +#process.MessageLogger = cms.Service("MessageLogger", +# destinations = cms.untracked.vstring("debugVH_tilted"), +# debugModules = cms.untracked.vstring("*"), +# categories = cms.untracked.vstring("VectorHitBuilderEDProducer","VectorHitBuilderAlgorithm","VectorHitsBuilderValidation","VectorHitBuilder"), +# debugVH_tilted = cms.untracked.PSet(threshold = cms.untracked.string("DEBUG"), +# DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(0)), +# default = cms.untracked.PSet(limit = cms.untracked.int32(0)), +# VectorHitBuilder = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# VectorHitBuilderEDProducer = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# VectorHitBuilderAlgorithm = cms.untracked.PSet(limit = cms.untracked.int32(-1)), +# VectorHitsBuilderValidation = cms.untracked.PSet(limit = cms.untracked.int32(-1)) +# ) +# ) + +# Analyzer +process.analysis = cms.EDAnalyzer('VectorHitsBuilderValidation', + src = cms.string("siPhase2Clusters"), + VH_acc = cms.InputTag("siPhase2VectorHits", "accepted"), + VH_rej = cms.InputTag("siPhase2VectorHits", "rejected"), + CPE = cms.ESInputTag("phase2StripCPEESProducer", "Phase2StripCPE"), + links = cms.InputTag("simSiPixelDigis", "Tracker"), + trackingParticleSrc = cms.InputTag('mix', 'MergedTrackTruth'), +) +process.TFileService = cms.Service('TFileService', + fileName = cms.string('file:VHs_validation.root') +) + + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.trackerlocalreco_step = cms.Path(process.trackerlocalreco+process.siPhase2VectorHits) +process.analysis_step = cms.Path(process.analysis) +process.RECOSIMoutput_step = cms.EndPath(process.RECOSIMoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.trackerlocalreco_step,process.RECOSIMoutput_step, process.analysis_step) + diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.cc b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.cc new file mode 100644 index 0000000000000..aac23a7052c51 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.cc @@ -0,0 +1,1137 @@ +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" +#include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.h" +#include "Geometry/CommonDetUnit/interface/StackGeomDet.h" +#include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "RecoLocalTracker/Records/interface/TkPhase2OTCPERecord.h" +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" + +VectorHitsBuilderValidation::VectorHitsBuilderValidation(const edm::ParameterSet& conf) + : cpeTag_(conf.getParameter("CPE")) { + srcClu_ = + consumes >(edm::InputTag(conf.getParameter("src"))); + VHacc_ = consumes(edm::InputTag(conf.getParameter("VH_acc"))); + VHrej_ = consumes(edm::InputTag(conf.getParameter("VH_rej"))); + siphase2OTSimLinksToken_ = consumes >(conf.getParameter("links")); + simHitsToken_ = consumes(edm::InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof")); + simTracksToken_ = consumes(edm::InputTag("g4SimHits")); + simVerticesToken_ = consumes(edm::InputTag("g4SimHits")); + trackingParticleToken_ = + consumes(conf.getParameter("trackingParticleSrc")); +} + +VectorHitsBuilderValidation::~VectorHitsBuilderValidation() {} + +void VectorHitsBuilderValidation::beginJob() { + edm::Service fs; + fs->file().cd("/"); + TFileDirectory td = fs->mkdir("Common"); + + //Create common ntuple + tree_ = td.make("VectorHits", "VectorHits"); + + // Create common graphs + TFileDirectory tdGloPos = td.mkdir("GlobalPositions"); + trackerLayoutRZ_[0] = tdGloPos.make(); + trackerLayoutRZ_[0]->SetName("RVsZ_Mixed"); + trackerLayoutRZ_[1] = tdGloPos.make(); + trackerLayoutRZ_[1]->SetName("RVsZ_Pixel"); + trackerLayoutRZ_[2] = tdGloPos.make(); + trackerLayoutRZ_[2]->SetName("RVsZ_Strip"); + trackerLayoutXY_[0] = tdGloPos.make(); + trackerLayoutXY_[0]->SetName("YVsX_Mixed"); + trackerLayoutXY_[1] = tdGloPos.make(); + trackerLayoutXY_[1]->SetName("YVsX_Pixel"); + trackerLayoutXY_[2] = tdGloPos.make(); + trackerLayoutXY_[2]->SetName("YVsX_Strip"); + trackerLayoutXYBar_ = tdGloPos.make(); + trackerLayoutXYBar_->SetName("YVsXBar"); + trackerLayoutXYEC_ = tdGloPos.make(); + trackerLayoutXYEC_->SetName("YVsXEC"); + + TFileDirectory tdLocPos = td.mkdir("LocalPositions"); + localPosXvsDeltaX_[0] = tdLocPos.make(); + localPosXvsDeltaX_[0]->SetName("localPosXvsDeltaX_Mixed"); + localPosXvsDeltaX_[1] = tdLocPos.make(); + localPosXvsDeltaX_[1]->SetName("localPosXvsDeltaX_Pixel"); + localPosXvsDeltaX_[2] = tdLocPos.make(); + localPosXvsDeltaX_[2]->SetName("localPosXvsDeltaX_Strip"); + localPosYvsDeltaY_[0] = tdLocPos.make(); + localPosYvsDeltaY_[0]->SetName("localPosYvsDeltaY_Mixed"); + localPosYvsDeltaY_[1] = tdLocPos.make(); + localPosYvsDeltaY_[1]->SetName("localPosYvsDeltaY_Pixel"); + localPosYvsDeltaY_[2] = tdLocPos.make(); + localPosYvsDeltaY_[2]->SetName("localPosYvsDeltaY_Strip"); + + //drawing VHs arrows + TFileDirectory tdArr = td.mkdir("Directions"); + + TFileDirectory tdWid = td.mkdir("CombinatorialStudies"); + ParallaxCorrectionRZ_ = + tdWid.make("ParallaxCorrectionRZ", "ParallaxCorrectionRZ", 100, 0., 300., 100., 0., 120.); + ParallaxCorrectionRZ_->SetName("ParallaxCorrectionFactor"); + VHaccLayer_ = tdWid.make("VHacceptedLayer", "VHacceptedLayer", 250, 0., 250.); + VHaccLayer_->SetName("VHaccepted"); + VHrejLayer_ = tdWid.make("VHrejectedLayer", "VHrejectedLayer", 250, 0., 250.); + VHrejLayer_->SetName("VHrejected"); + VHaccTrueLayer_ = tdWid.make("VHaccTrueLayer", "VHaccTrueLayer", 250, 0., 250.); + VHaccTrueLayer_->SetName("VHaccepted_true"); + VHrejTrueLayer_ = tdWid.make("VHrejTrueLayer", "VHrejTrueLayer", 250, 0., 250.); + VHrejTrueLayer_->SetName("VHrejected_true"); + VHaccTrue_signal_Layer_ = tdWid.make("VHaccTrueSignalLayer", "VHaccTrueSignalLayer", 250, 0., 250.); + VHaccTrue_signal_Layer_->SetName("VHaccepted_true_signal"); + VHrejTrue_signal_Layer_ = tdWid.make("VHrejTrueSignalLayer", "VHrejTrueSignalLayer", 250, 0., 250.); + VHrejTrue_signal_Layer_->SetName("VHrejected_true_signal"); +} + +void VectorHitsBuilderValidation::endJob() {} + +void VectorHitsBuilderValidation::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) { + // Get the needed objects + + // Get the clusters + edm::Handle clusters; + event.getByToken(srcClu_, clusters); + + // Get the vector hits + edm::Handle vhsAcc; + event.getByToken(VHacc_, vhsAcc); + + edm::Handle vhsRej; + event.getByToken(VHrej_, vhsRej); + + // load the cpe via the eventsetup + edm::ESHandle > cpeHandle; + eventSetup.get().get(cpeTag_, cpeHandle); + cpe_ = cpeHandle.product(); + + // Get the Phase2 DigiSimLink + edm::Handle > siphase2SimLinks; + event.getByToken(siphase2OTSimLinksToken_, siphase2SimLinks); + + // Get the SimHits + edm::Handle simHitsRaw; + event.getByToken(simHitsToken_, simHitsRaw); + + // Get the SimTracks + edm::Handle simTracksRaw; + event.getByToken(simTracksToken_, simTracksRaw); + + // Get the SimVertex + edm::Handle simVertices; + event.getByToken(simVerticesToken_, simVertices); + + // Get the geometry + edm::ESHandle geomHandle; + eventSetup.get().get(geomHandle); + tkGeom_ = &(*geomHandle); + + // Get the Topology + edm::ESHandle tTopoHandle; + eventSetup.get().get(tTopoHandle); + tkTopo_ = tTopoHandle.product(); + + edm::ESHandle magFieldHandle; + eventSetup.get().get(magFieldHandle); + magField_ = magFieldHandle.product(); + + //Tracking Particle collection + edm::Handle TPCollectionH; + event.getByToken(trackingParticleToken_, TPCollectionH); + + auto clusterTPList = std::make_unique(TPCollectionH); + std::map, TrackingParticleRef> mapping; + + for (TrackingParticleCollection::size_type itp = 0; itp < TPCollectionH.product()->size(); ++itp) { + TrackingParticleRef trackingParticle(TPCollectionH, itp); + EncodedEventId eid(trackingParticle->eventId()); + for (std::vector::const_iterator itrk = trackingParticle->g4Track_begin(); + itrk != trackingParticle->g4Track_end(); + ++itrk) { + std::pair trkid(itrk->trackId(), eid); + LogTrace("VectorHitsBuilderValidation") + << "creating map for id: " << trkid.first << " with tp: " << trackingParticle.key(); + mapping.insert(std::make_pair(trkid, trackingParticle)); + } + } + + //set up for tree + int eventNum; + int layer; + int module_id; + int module_number; + int module_type; //1: pixel, 2: strip + int VHacc = 0.0; + int VHrej = 0.0; + int vh_isTrue; + + float x_global, y_global, z_global; + float vh_x_local, vh_y_local; + float vh_x_le, vh_y_le; + float curvature, phi; + float QOverPT, QOverP; + float chi2; + + int low_tp_id, upp_tp_id; + float vh_sim_trackPt; + float sim_x_local, sim_y_local; + float sim_x_global, sim_y_global, sim_z_global; + float low_x_global, low_y_global, low_z_global; + float upp_x_global, upp_y_global, upp_z_global; + float low_xx_global_err, low_yy_global_err, low_zz_global_err; + float low_xy_global_err, low_zx_global_err, low_zy_global_err; + float upp_xx_global_err, upp_yy_global_err, upp_zz_global_err; + float upp_xy_global_err, upp_zx_global_err, upp_zy_global_err; + float deltaXVHSimHits, deltaYVHSimHits; + int multiplicity; + float width, deltaXlocal; + unsigned int processType(99); + + tree_->Branch("event", &eventNum, "eventNum/I"); + tree_->Branch("accepted", &VHacc, "VHacc/I"); + tree_->Branch("rejected", &VHrej, "VHrej/I"); + tree_->Branch("layer", &layer, "layer/I"); + tree_->Branch("module_id", &module_id, "module_id/I"); + tree_->Branch("module_type", &module_type, "module_type/I"); + tree_->Branch("module_number", &module_number, "module_number/I"); + tree_->Branch("vh_isTrue", &vh_isTrue, "vh_isTrue/I"); + tree_->Branch("x_global", &x_global, "x_global/F"); + tree_->Branch("y_global", &y_global, "y_global/F"); + tree_->Branch("z_global", &z_global, "z_global/F"); + tree_->Branch("vh_x_local", &vh_x_local, "vh_x_local/F"); + tree_->Branch("vh_y_local", &vh_y_local, "vh_y_local/F"); + tree_->Branch("vh_x_lError", &vh_x_le, "vh_x_le/F"); + tree_->Branch("vh_y_lError", &vh_y_le, "vh_y_le/F"); + tree_->Branch("curvature", &curvature, "curvature/F"); + tree_->Branch("chi2", &chi2, "chi2/F"); + tree_->Branch("phi", &phi, "phi/F"); + tree_->Branch("QOverP", &QOverP, "QOverP/F"); + tree_->Branch("QOverPT", &QOverPT, "QOverPT/F"); + tree_->Branch("low_tp_id", &low_tp_id, "low_tp_id/I"); + tree_->Branch("upp_tp_id", &upp_tp_id, "upp_tp_id/I"); + tree_->Branch("vh_sim_trackPt", &vh_sim_trackPt, "vh_sim_trackPt/F"); + tree_->Branch("sim_x_local", &sim_x_local, "sim_x_local/F"); + tree_->Branch("sim_y_local", &sim_y_local, "sim_y_local/F"); + tree_->Branch("sim_x_global", &sim_x_global, "sim_x_global/F"); + tree_->Branch("sim_y_global", &sim_y_global, "sim_y_global/F"); + tree_->Branch("sim_z_global", &sim_z_global, "sim_z_global/F"); + tree_->Branch("low_x_global", &low_x_global, "low_x_global/F"); + tree_->Branch("low_y_global", &low_y_global, "low_y_global/F"); + tree_->Branch("low_z_global", &low_z_global, "low_z_global/F"); + tree_->Branch("low_xx_global_err", &low_xx_global_err, "low_xx_global_err/F"); + tree_->Branch("low_yy_global_err", &low_yy_global_err, "low_yy_global_err/F"); + tree_->Branch("low_zz_global_err", &low_zz_global_err, "low_zz_global_err/F"); + tree_->Branch("low_xy_global_err", &low_xy_global_err, "low_xy_global_err/F"); + tree_->Branch("low_zx_global_err", &low_zx_global_err, "low_zx_global_err/F"); + tree_->Branch("low_zy_global_err", &low_zy_global_err, "low_zy_global_err/F"); + tree_->Branch("upp_x_global", &upp_x_global, "upp_x_global/F"); + tree_->Branch("upp_y_global", &upp_y_global, "upp_y_global/F"); + tree_->Branch("upp_z_global", &upp_z_global, "upp_z_global/F"); + tree_->Branch("upp_xx_global_err", &upp_xx_global_err, "upp_xx_global_err/F"); + tree_->Branch("upp_yy_global_err", &upp_yy_global_err, "upp_yy_global_err/F"); + tree_->Branch("upp_zz_global_err", &upp_zz_global_err, "upp_zz_global_err/F"); + tree_->Branch("upp_xy_global_err", &upp_xy_global_err, "upp_xy_global_err/F"); + tree_->Branch("upp_zx_global_err", &upp_zx_global_err, "upp_zx_global_err/F"); + tree_->Branch("upp_zy_global_err", &upp_zy_global_err, "upp_zy_global_err/F"); + tree_->Branch("deltaXVHSimHits", &deltaXVHSimHits, "deltaXVHSimHits/F"); + tree_->Branch("deltaYVHSimHits", &deltaYVHSimHits, "deltaYVHSimHits/F"); + tree_->Branch("multiplicity", &multiplicity, "multiplicity/I"); + tree_->Branch("width", &width, "width/F"); + tree_->Branch("deltaXlocal", &deltaXlocal, "deltaXlocal/F"); + tree_->Branch("processType", &processType, "processType/i"); + + // Rearrange the simTracks for ease of use + SimTracksMap simTracks; + for (const auto& simTrackIt : *simTracksRaw.product()) + simTracks.emplace(std::pair(simTrackIt.trackId(), simTrackIt)); + + // Rearrange the simHits by detUnit for ease of use + SimHitsMap simHitsDetUnit; + SimHitsMap simHitsTrackId; + for (const auto& simHitIt : *simHitsRaw.product()) { + SimHitsMap::iterator simHitsDetUnitIt(simHitsDetUnit.find(simHitIt.detUnitId())); + if (simHitsDetUnitIt == simHitsDetUnit.end()) { + std::pair newIt(simHitsDetUnit.insert( + std::pair >(simHitIt.detUnitId(), std::vector()))); + simHitsDetUnitIt = newIt.first; + } + simHitsDetUnitIt->second.push_back(simHitIt); + + SimHitsMap::iterator simHitsTrackIdIt(simHitsTrackId.find(simHitIt.trackId())); + if (simHitsTrackIdIt == simHitsTrackId.end()) { + std::pair newIt(simHitsTrackId.insert( + std::pair >(simHitIt.trackId(), std::vector()))); + simHitsTrackIdIt = newIt.first; + } + simHitsTrackIdIt->second.push_back(simHitIt); + } + + //Printout outer tracker clusters in the event + for (const auto& DSViter : *clusters) { + unsigned int rawid(DSViter.detId()); + DetId detId(rawid); + const GeomDetUnit* geomDetUnit(tkGeom_->idToDetUnit(detId)); + const PixelGeomDetUnit* theGeomDet = dynamic_cast(geomDetUnit); + for (const auto& clustIt : DSViter) { + auto&& lparams = cpe_->localParameters(clustIt, *theGeomDet); + Global3DPoint gparams = theGeomDet->surface().toGlobal(lparams.first); + LogTrace("VectorHitsBuilderValidation") << "phase2 OT clusters: " << gparams << " DetId: " << rawid; + } + } + + for (const auto& DSViter : *vhsAcc) { + for (const auto& vhIt : DSViter) { + LogTrace("VectorHitsBuilderValidation") << "accepted VH: " << vhIt; + } + } + for (const auto& DSViter : *vhsRej) { + for (const auto& vhIt : DSViter) { + LogTrace("VectorHitsBuilderValidation") << "rejected VH: " << vhIt; + } + } + // Validation + eventNum = event.id().event(); + + unsigned int nVHsTot(0), nVHsPSTot(0), nVHs2STot(0); + std::vector glVHs; + std::vector dirVHs; + std::vector detIds; + + // Loop over modules + for (const auto DSViter : *vhsAcc) { + // Get the detector unit's id + unsigned int rawid(DSViter.detId()); + module_id = rawid; + DetId detId(rawid); + + module_number = getModuleNumber(detId); + layer = getLayerNumber(detId); + + LogDebug("VectorHitsBuilderValidation") << "Layer: " << layer << " det id" << rawid << std::endl; + + // Get the geometry of the tracker + const GeomDet* geomDet(tkGeom_->idToDet(detId)); + if (!geomDet) + break; + + // Create histograms for the layer if they do not yet exist + std::map::iterator histogramLayer(histograms_.find(layer)); + if (histogramLayer == histograms_.end()) + histogramLayer = createLayerHistograms(layer); + // Number of clusters + unsigned int nVHsPS(0), nVHs2S(0); + + LogDebug("VectorHitsBuilderValidation") << "DSViter size: " << DSViter.size(); + + // Loop over the vhs in the detector unit + for (const auto& vhIt : DSViter) { + // vh variables + if (vhIt.isValid()) { + LogDebug("VectorHitsBuilderValidation") << " vh analyzing ..."; + chi2 = vhIt.chi2(); + LogTrace("VectorHitsBuilderValidation") << "VH chi2 " << chi2 << std::endl; + + Local3DPoint localPosVH = vhIt.localPosition(); + vh_x_local = localPosVH.x(); + vh_y_local = localPosVH.y(); + LogTrace("VectorHitsBuilderValidation") << "local VH position " << localPosVH << std::endl; + + LocalError localErrVH = vhIt.localPositionError(); + vh_x_le = localErrVH.xx(); + vh_y_le = localErrVH.yy(); + LogTrace("VectorHitsBuilderValidation") << "local VH error " << localErrVH << std::endl; + + Global3DPoint globalPosVH = geomDet->surface().toGlobal(localPosVH); + x_global = globalPosVH.x(); + y_global = globalPosVH.y(); + z_global = globalPosVH.z(); + glVHs.push_back(globalPosVH); + LogTrace("VectorHitsBuilderValidation") << " global VH position " << globalPosVH << std::endl; + + Local3DVector localDirVH = vhIt.localDirection(); + LogTrace("VectorHitsBuilderValidation") << "local VH direction " << localDirVH << std::endl; + + VectorHit vh = vhIt; + Global3DVector globalDirVH = vh.globalDirectionVH(); + dirVHs.push_back(globalDirVH); + LogTrace("VectorHitsBuilderValidation") << "global VH direction " << globalDirVH << std::endl; + + // Fill the position histograms + trackerLayoutRZ_[0]->SetPoint(nVHsTot, globalPosVH.z(), globalPosVH.perp()); + trackerLayoutXY_[0]->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y()); + + if (layer < 100) + trackerLayoutXYBar_->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y()); + else + trackerLayoutXYEC_->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y()); + + histogramLayer->second.localPosXY[0]->SetPoint(nVHsTot, vh_x_local, vh_y_local); + histogramLayer->second.globalPosXY[0]->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y()); + + localPosXvsDeltaX_[0]->SetPoint(nVHsTot, vh_x_local, localDirVH.x()); + localPosYvsDeltaY_[0]->SetPoint(nVHsTot, vh_y_local, localDirVH.y()); + + // Pixel module + const StackGeomDet* stackDet = dynamic_cast(geomDet); + const PixelGeomDetUnit* geomDetLower = dynamic_cast(stackDet->lowerDet()); + DetId lowerDetId = stackDet->lowerDet()->geographicalId(); + DetId upperDetId = stackDet->upperDet()->geographicalId(); + + TrackerGeometry::ModuleType mType = tkGeom_->getDetectorType(lowerDetId); + module_type = 0; + if (mType == TrackerGeometry::ModuleType::Ph2PSP) { + module_type = 1; + trackerLayoutRZ_[1]->SetPoint(nVHsPSTot, globalPosVH.z(), globalPosVH.perp()); + trackerLayoutXY_[1]->SetPoint(nVHsPSTot, globalPosVH.x(), globalPosVH.y()); + + histogramLayer->second.localPosXY[1]->SetPoint(nVHsPSTot, vh_x_local, vh_y_local); + histogramLayer->second.globalPosXY[1]->SetPoint(nVHsPSTot, globalPosVH.x(), globalPosVH.y()); + + localPosXvsDeltaX_[1]->SetPoint(nVHsPSTot, vh_x_local, localDirVH.x()); + localPosYvsDeltaY_[1]->SetPoint(nVHsPSTot, vh_y_local, localDirVH.y()); + + ++nVHsPS; + ++nVHsPSTot; + } + + // Strip module + else if (mType == TrackerGeometry::ModuleType::Ph2SS) { + module_type = 2; + trackerLayoutRZ_[2]->SetPoint(nVHs2STot, globalPosVH.z(), globalPosVH.perp()); + trackerLayoutXY_[2]->SetPoint(nVHs2STot, globalPosVH.x(), globalPosVH.y()); + + histogramLayer->second.localPosXY[2]->SetPoint(nVHs2STot, vh_x_local, vh_y_local); + histogramLayer->second.globalPosXY[2]->SetPoint(nVHs2STot, globalPosVH.x(), globalPosVH.y()); + + localPosXvsDeltaX_[2]->SetPoint(nVHs2STot, vh_x_local, localDirVH.x()); + localPosYvsDeltaY_[2]->SetPoint(nVHs2STot, vh_y_local, localDirVH.y()); + + ++nVHs2S; + ++nVHs2STot; + } else if (mType == TrackerGeometry::ModuleType::Ph2PSS) { + edm::LogError("VectorHitsBuilderValidation") << "module type " << module_type << " should never happen!"; + } + LogTrace("VectorHitsBuilderValidation") << "module type " << module_type << std::endl; + + // get the geomDetUnit of the clusters + low_x_global = vhIt.lowerGlobalPos().x(); + low_y_global = vhIt.lowerGlobalPos().y(); + low_z_global = vhIt.lowerGlobalPos().z(); + upp_x_global = vhIt.upperGlobalPos().x(); + upp_y_global = vhIt.upperGlobalPos().y(); + upp_z_global = vhIt.upperGlobalPos().z(); + + low_xx_global_err = vhIt.lowerGlobalPosErr().cxx(); + low_yy_global_err = vhIt.lowerGlobalPosErr().cyy(); + low_zz_global_err = vhIt.lowerGlobalPosErr().czz(); + low_xy_global_err = vhIt.lowerGlobalPosErr().cyx(); + low_zx_global_err = vhIt.lowerGlobalPosErr().czx(); + low_zy_global_err = vhIt.lowerGlobalPosErr().czy(); + + upp_xx_global_err = vhIt.upperGlobalPosErr().cxx(); + upp_yy_global_err = vhIt.upperGlobalPosErr().cyy(); + upp_zz_global_err = vhIt.upperGlobalPosErr().czz(); + upp_xy_global_err = vhIt.upperGlobalPosErr().cyx(); + upp_zx_global_err = vhIt.upperGlobalPosErr().czx(); + upp_zy_global_err = vhIt.upperGlobalPosErr().czy(); + + LogDebug("VectorHitsBuilderValidation") << "print Clusters into the VH:" << std::endl; + printCluster(geomDetLower, vhIt.lowerClusterRef()); + LogTrace("VectorHitsBuilderValidation") << "\t global pos lower " << vhIt.lowerGlobalPos() << std::endl; + LogTrace("VectorHitsBuilderValidation") + << "\t global posErr lower " << vhIt.lowerGlobalPosErr().cxx() << std::endl; + const GeomDetUnit* geomDetUpper = stackDet->upperDet(); + printCluster(geomDetUpper, vhIt.upperClusterRef()); + LogTrace("VectorHitsBuilderValidation") << "\t global pos upper " << vhIt.upperGlobalPos() << std::endl; + + //comparison with SIM hits + LogDebug("VectorHitsBuilderValidation") << "comparison Clusters with sim hits ... " << std::endl; + std::vector clusterSimTrackIds; + std::vector clusterSimTrackIdsUpp; + std::set > simTkIds; + const GeomDetUnit* geomDetUnit_low(tkGeom_->idToDetUnit(lowerDetId)); + LogTrace("VectorHitsBuilderValidation") << " lowerDetID : " << lowerDetId.rawId(); + const GeomDetUnit* geomDetUnit_upp(tkGeom_->idToDetUnit(upperDetId)); + LogTrace("VectorHitsBuilderValidation") << " upperDetID : " << upperDetId.rawId(); + + for (unsigned int istr(0); istr < (*(vhIt.lowerClusterRef().cluster_phase2OT())).size(); ++istr) { + uint32_t channel = + Phase2TrackerDigi::pixelToChannel((*(vhIt.lowerClusterRef().cluster_phase2OT())).firstRow() + istr, + (*(vhIt.lowerClusterRef().cluster_phase2OT())).column()); + unsigned int LowerSimTrackId(getSimTrackId(siphase2SimLinks, lowerDetId, channel)); + std::vector > trkid( + getSimTrackIds(siphase2SimLinks, lowerDetId, channel)); + if (trkid.size() == 0) + continue; + clusterSimTrackIds.push_back(LowerSimTrackId); + simTkIds.insert(trkid.begin(), trkid.end()); + LogTrace("VectorHitsBuilderValidation") << "LowerSimTrackId " << LowerSimTrackId << std::endl; + } + // In the case of PU, we need the TPs to find the proper SimTrackID + for (const auto& iset : simTkIds) { + auto ipos = mapping.find(iset); + if (ipos != mapping.end()) { + LogTrace("VectorHitsBuilderValidation") << "lower cluster in detid: " << lowerDetId.rawId() + << " from tp: " << ipos->second.key() << " " << iset.first; + LogTrace("VectorHitsBuilderValidation") << "with pt(): " << (*ipos->second).pt(); + low_tp_id = ipos->second.key(); + vh_sim_trackPt = (*ipos->second).pt(); + } + } + + simTkIds.clear(); + for (unsigned int istr(0); istr < (*(vhIt.upperClusterRef().cluster_phase2OT())).size(); ++istr) { + uint32_t channel = + Phase2TrackerDigi::pixelToChannel((*(vhIt.upperClusterRef().cluster_phase2OT())).firstRow() + istr, + (*(vhIt.upperClusterRef().cluster_phase2OT())).column()); + unsigned int UpperSimTrackId(getSimTrackId(siphase2SimLinks, upperDetId, channel)); + std::vector > trkid( + getSimTrackIds(siphase2SimLinks, upperDetId, channel)); + if (trkid.size() == 0) + continue; + clusterSimTrackIdsUpp.push_back(UpperSimTrackId); + simTkIds.insert(trkid.begin(), trkid.end()); + LogTrace("VectorHitsBuilderValidation") << "UpperSimTrackId " << UpperSimTrackId << std::endl; + } + // In the case of PU, we need the TPs to find the proper SimTrackID + for (const auto& iset : simTkIds) { + auto ipos = mapping.find(iset); + if (ipos != mapping.end()) { + LogTrace("VectorHitsBuilderValidation") + << "upper cluster in detid: " << upperDetId.rawId() << " from tp: " << ipos->second.key() << " " + << iset.first << std::endl; + upp_tp_id = ipos->second.key(); + } + } + //compute if the vhits is 'true' or 'false' and save sim pT + std::pair istrue = isTrue(vhIt, siphase2SimLinks, detId); + vh_isTrue = 0; + if (istrue.first) { + vh_isTrue = 1; + } + + // loop over all simHits + unsigned int totalSimHits(0); + unsigned int primarySimHits(0); + unsigned int otherSimHits(0); + + for (const auto& hitIt : *simHitsRaw) { + if (hitIt.detUnitId() == geomDetLower->geographicalId()) { + //check clusters track id compatibility + if (std::find(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), hitIt.trackId()) != + clusterSimTrackIds.end()) { + Local3DPoint localPosHit(hitIt.localPosition()); + sim_x_local = localPosHit.x(); + sim_y_local = localPosHit.y(); + + deltaXVHSimHits = vh_x_local - sim_x_local; + deltaYVHSimHits = vh_y_local - sim_y_local; + + Global3DPoint globalPosHit = geomDetLower->surface().toGlobal(localPosHit); + sim_x_global = globalPosHit.x(); + sim_y_global = globalPosHit.y(); + sim_z_global = globalPosHit.z(); + + histogramLayer->second.deltaXVHSimHits[0]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits[0]->Fill(vh_y_local - sim_y_local); + + // Pixel module + if (layer == 1 || layer == 2 || layer == 3) { + histogramLayer->second.deltaXVHSimHits[1]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits[1]->Fill(vh_y_local - sim_y_local); + } + // Strip module + else if (layer == 4 || layer == 5 || layer == 6) { + histogramLayer->second.deltaXVHSimHits[2]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits[2]->Fill(vh_y_local - sim_y_local); + } + + ++totalSimHits; + + std::map::const_iterator simTrackIt(simTracks.find(hitIt.trackId())); + if (simTrackIt == simTracks.end()) + continue; + + // Primary particles only + processType = hitIt.processType(); + + if (simTrackIt->second.vertIndex() == 0 and + (processType == 2 || processType == 7 || processType == 9 || processType == 11 || processType == 13 || + processType == 15)) { + histogramLayer->second.deltaXVHSimHits_P[0]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits_P[0]->Fill(vh_y_local - sim_y_local); + + // Pixel module + if (layer == 1 || layer == 2 || layer == 3) { + histogramLayer->second.deltaXVHSimHits_P[1]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits_P[1]->Fill(vh_y_local - sim_y_local); + } + // Strip module + else if (layer == 4 || layer == 5 || layer == 6) { + histogramLayer->second.deltaXVHSimHits_P[2]->Fill(vh_x_local - sim_x_local); + histogramLayer->second.deltaYVHSimHits_P[2]->Fill(vh_y_local - sim_y_local); + } + + ++primarySimHits; + } + + otherSimHits = totalSimHits - primarySimHits; + + histogramLayer->second.totalSimHits->Fill(totalSimHits); + histogramLayer->second.primarySimHits->Fill(primarySimHits); + histogramLayer->second.otherSimHits->Fill(otherSimHits); + } + } + } // loop simhits + + nVHsTot++; + + //****************************** + //combinatorial studies : not filling if more than 1 VH has been produced + //****************************** + multiplicity = DSViter.size(); + if (DSViter.size() > 1) { + LogTrace("VectorHitsBuilderValidation") << " not filling if more than 1 VH has been produced"; + width = -100; + deltaXlocal = -100; + tree_->Fill(); + continue; + } + + //curvature + GlobalPoint center(0.0, 0.0, 0.0); + curvature = vh.curvature(); + phi = vh.phi(); + QOverPT = vh.transverseMomentum(magField_->inTesla(center).z()); + QOverP = vh.momentum(magField_->inTesla(center).z()); + histogramLayer->second.curvature->Fill(curvature); + + //stub width + + auto&& lparamsUpp = cpe_->localParameters(*vhIt.upperClusterRef().cluster_phase2OT(), *geomDetUnit_upp); + LogTrace("VectorHitsBuilderValidation") << " upper local pos (in its system of reference):" << lparamsUpp.first; + Global3DPoint gparamsUpp = geomDetUnit_upp->surface().toGlobal(lparamsUpp.first); + LogTrace("VectorHitsBuilderValidation") << " upper global pos :" << gparamsUpp; + Local3DPoint lparamsUppInLow = geomDetUnit_low->surface().toLocal(gparamsUpp); + LogTrace("VectorHitsBuilderValidation") << " upper local pos (in low system of reference):" << lparamsUppInLow; + auto&& lparamsLow = cpe_->localParameters(*vhIt.lowerClusterRef().cluster_phase2OT(), *geomDetUnit_low); + LogTrace("VectorHitsBuilderValidation") << " lower local pos (in its system of reference):" << lparamsLow.first; + Global3DPoint gparamsLow = geomDetUnit_low->surface().toGlobal(lparamsLow.first); + LogTrace("VectorHitsBuilderValidation") << " lower global pos :" << gparamsLow; + + deltaXlocal = lparamsUppInLow.x() - lparamsLow.first.x(); + histogramLayer->second.deltaXlocal->Fill(deltaXlocal); + LogTrace("VectorHitsBuilderValidation") << " deltaXlocal : " << deltaXlocal; + + double parallCorr = 0.0; + + Global3DPoint origin(0, 0, 0); + GlobalVector gV = gparamsLow - origin; + LocalVector lV = geomDetUnit_low->surface().toLocal(gV); + LocalVector lV_norm = lV / lV.z(); + parallCorr = lV_norm.x() * lparamsUppInLow.z(); + LogTrace("VectorHitsBuilderValidation") << " parallalex correction:" << parallCorr; + + double lpos_upp_corr = 0.0; + double lpos_low_corr = 0.0; + if (lparamsUpp.first.x() > lparamsLow.first.x()) { + if (lparamsUpp.first.x() > 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = lparamsUpp.first.x() - fabs(parallCorr); + } + if (lparamsUpp.first.x() < 0) { + lpos_low_corr = lparamsLow.first.x() + fabs(parallCorr); + lpos_upp_corr = lparamsUpp.first.x(); + } + } else if (lparamsUpp.first.x() < lparamsLow.first.x()) { + if (lparamsUpp.first.x() > 0) { + lpos_low_corr = lparamsLow.first.x() - fabs(parallCorr); + lpos_upp_corr = lparamsUpp.first.x(); + } + if (lparamsUpp.first.x() < 0) { + lpos_low_corr = lparamsLow.first.x(); + lpos_upp_corr = lparamsUpp.first.x() + fabs(parallCorr); + } + } else { + if (lparamsUpp.first.x() > 0) { + lpos_upp_corr = lparamsUpp.first.x() - fabs(parallCorr); + lpos_low_corr = lparamsLow.first.x(); + } + if (lparamsUpp.first.x() < 0) { + lpos_upp_corr = lparamsUpp.first.x() + fabs(parallCorr); + lpos_low_corr = lparamsLow.first.x(); + } + } + + LogDebug("VectorHitsBuilderValidation") << " \t local pos upper corrected (x):" << lpos_upp_corr << std::endl; + LogDebug("VectorHitsBuilderValidation") << " \t local pos lower corrected (x):" << lpos_low_corr << std::endl; + + width = lpos_low_corr - lpos_upp_corr; + histogramLayer->second.width->Fill(width); + LogTrace("VectorHitsBuilderValidation") << " width:" << width; + + tree_->Fill(); + + } // vh valid + + } // loop vhs + + if (nVHsPS) + histogramLayer->second.numberVHsPS->Fill(nVHsPS); + if (nVHs2S) + histogramLayer->second.numberVHs2S->Fill(nVHs2S); + LogTrace("VectorHitsBuilderValidation") + << "nVHsPS for this layer : " << nVHsPS << ", nVHs2S for this layer : " << nVHs2S << std::endl; + } + + CreateVHsXYGraph(glVHs, dirVHs); + CreateVHsRZGraph(glVHs, dirVHs); + + int VHaccTrue = 0.0; + int VHaccFalse = 0.0; + int VHrejTrue = 0.0; + int VHrejFalse = 0.0; + int VHaccTrue_signal = 0.0; + int VHrejTrue_signal = 0.0; + + // Loop over modules + for (const auto& DSViter : *vhsAcc) { + unsigned int rawid(DSViter.detId()); + DetId detId(rawid); + int layerAcc = getLayerNumber(detId); + LogTrace("VectorHitsBuilderValidation") << "acc Layer: " << layerAcc << " det id" << rawid << std::endl; + for (const auto& vhIt : DSViter) { + if (vhIt.isValid()) { + VHaccLayer_->Fill(layerAcc); + VHacc++; + + //compute if the vhits is 'true' or 'false' + std::pair istrue = isTrue(vhIt, siphase2SimLinks, detId); + if (istrue.first) { + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'true' vhit."; + VHaccTrueLayer_->Fill(layerAcc); + VHaccTrue++; + + //saving info of 'signal' track + std::map::const_iterator simTrackIt(simTracks.find(istrue.second)); + if (simTrackIt == simTracks.end()) + continue; + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is associated with SimTrackId: " << istrue.second; + LogTrace("VectorHitsBuilderValidation") << "the SimTrack has pt: " << simTrackIt->second.momentum().pt(); + if (simTrackIt->second.momentum().pt() > 1) { + VHaccTrue_signal_Layer_->Fill(layerAcc); + LogTrace("VectorHitsBuilderValidation") << "the vectorhit belongs to signal"; + VHaccTrue_signal++; + } + + } else { + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'false' vhit."; + VHaccFalse++; + } + } + } + } + + for (const auto& DSViter : *vhsRej) { + unsigned int rawid(DSViter.detId()); + DetId detId(rawid); + int layerRej = getLayerNumber(detId); + LogTrace("VectorHitsBuilderValidation") << "rej Layer: " << layerRej << " det id" << rawid << std::endl; + for (const auto& vhIt : DSViter) { + VHrejLayer_->Fill(layerRej); + VHrej++; + + //compute if the vhits is 'true' or 'false' + std::pair istrue = isTrue(vhIt, siphase2SimLinks, detId); + if (istrue.first) { + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'true' vhit."; + VHrejTrueLayer_->Fill(layerRej); + VHrejTrue++; + + //saving info of 'signal' track + std::map::const_iterator simTrackIt(simTracks.find(istrue.second)); + if (simTrackIt == simTracks.end()) + continue; + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is associated with SimTrackId: " << istrue.second; + LogTrace("VectorHitsBuilderValidation") << "the SimTrack has pt: " << simTrackIt->second.momentum().pt(); + if (simTrackIt->second.momentum().pt() > 1) { + VHrejTrue_signal_Layer_->Fill(layerRej); + LogTrace("VectorHitsBuilderValidation") << "the vectorhit belongs to signal"; + VHrejTrue_signal++; + } + + } else { + LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'false' vhit."; + VHrejFalse++; + } + } + } + + int VHtot = VHacc + VHrej; + LogTrace("VectorHitsBuilderValidation") + << "VH total: " << VHtot << " with " << VHacc << " VHs accepted and " << VHrej << " VHs rejected."; + LogTrace("VectorHitsBuilderValidation") + << "of the VH accepted, there are " << VHaccTrue << " true and " << VHaccFalse << " false."; + LogTrace("VectorHitsBuilderValidation") + << "of the VH rejected, there are " << VHrejTrue << " true and " << VHrejFalse << " false."; + LogTrace("VectorHitsBuilderValidation") + << "of the true VH , there are " << VHaccTrue_signal << " accepted belonging to signal and " + << VHrejTrue_signal << " rejected belonging to signal."; + + // CreateWindowCorrGraph(); +} + +// Check if the vector hit is true (both clusters are formed from the same SimTrack +std::pair VectorHitsBuilderValidation::isTrue( + const VectorHit vh, const edm::Handle >& siphase2SimLinks, DetId& detId) const { + const GeomDet* geomDet(tkGeom_->idToDet(detId)); + const StackGeomDet* stackDet = dynamic_cast(geomDet); + const GeomDetUnit* geomDetLower = stackDet->lowerDet(); + const GeomDetUnit* geomDetUpper = stackDet->upperDet(); + + std::vector lowClusterSimTrackIds; + + for (unsigned int istr(0); istr < (*(vh.lowerClusterRef().cluster_phase2OT())).size(); ++istr) { + uint32_t channel = Phase2TrackerDigi::pixelToChannel((*(vh.lowerClusterRef().cluster_phase2OT())).firstRow() + istr, + (*(vh.lowerClusterRef().cluster_phase2OT())).column()); + DetId detIdCluster = geomDetLower->geographicalId(); + unsigned int simTrackId(getSimTrackId(siphase2SimLinks, detIdCluster, channel)); + LogTrace("VectorHitsBuilderValidation") << "LowerSimTrackId " << simTrackId << std::endl; + std::vector > trkid(getSimTrackIds(siphase2SimLinks, detIdCluster, channel)); + if (trkid.size() == 0) + continue; + lowClusterSimTrackIds.push_back(simTrackId); + } + + std::vector::iterator it_simTrackUpper; + + for (unsigned int istr(0); istr < (*(vh.upperClusterRef().cluster_phase2OT())).size(); ++istr) { + uint32_t channel = Phase2TrackerDigi::pixelToChannel((*(vh.upperClusterRef().cluster_phase2OT())).firstRow() + istr, + (*(vh.upperClusterRef().cluster_phase2OT())).column()); + DetId detIdCluster = geomDetUpper->geographicalId(); + unsigned int simTrackId(getSimTrackId(siphase2SimLinks, detIdCluster, channel)); + LogTrace("VectorHitsBuilderValidation") << "UpperSimTrackId " << simTrackId << std::endl; + std::vector > trkid(getSimTrackIds(siphase2SimLinks, detIdCluster, channel)); + if (trkid.size() == 0) + continue; + it_simTrackUpper = std::find(lowClusterSimTrackIds.begin(), lowClusterSimTrackIds.end(), simTrackId); + if (it_simTrackUpper != lowClusterSimTrackIds.end()) { + LogTrace("VectorHitsBuilderValidation") << " UpperSimTrackId found in lowClusterSimTrackIds "; + return std::make_pair(true, simTrackId); + } + } + return std::make_pair(false, 0); +} + +// Create the histograms +std::map::iterator VectorHitsBuilderValidation::createLayerHistograms(unsigned int ival) { + std::ostringstream fname1, fname2; + + edm::Service fs; + fs->file().cd("/"); + + std::string tag; + unsigned int id; + if (ival < 100) { + id = ival; + fname1 << "Barrel"; + fname2 << "Layer_" << id; + tag = "_layer_"; + } else { + int side = ival / 100; + id = ival - side * 100; + fname1 << "EndCap_Side_" << side; + fname2 << "Disc_" << id; + tag = "_disc_"; + } + + TFileDirectory td1 = fs->mkdir(fname1.str().c_str()); + TFileDirectory td = td1.mkdir(fname2.str().c_str()); + + VHHistos local_histos; + + std::ostringstream histoName; + + /* + * Number of clusters + */ + + histoName.str(""); + histoName << "Number_VHs_PS" << tag.c_str() << id; + local_histos.numberVHsPS = td.make(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.); + local_histos.numberVHsPS->SetFillColor(kAzure + 7); + + histoName.str(""); + histoName << "Number_VHs_2S" << tag.c_str() << id; + local_histos.numberVHs2S = td.make(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.); + local_histos.numberVHs2S->SetFillColor(kOrange - 3); + + histoName.str(""); + histoName << "Number_VHs_Mixed" << tag.c_str() << id; + local_histos.numberVHsMixed = td.make(histoName.str().c_str(), histoName.str().c_str()); + local_histos.numberVHsMixed->Add(local_histos.numberVHsPS); + local_histos.numberVHsMixed->Add(local_histos.numberVHs2S); + + /* + * Local and Global positions + */ + + histoName.str(""); + histoName << "Local_Position_XY_Mixed" << tag.c_str() << id; + local_histos.localPosXY[0] = td.make(); + local_histos.localPosXY[0]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Local_Position_XY_PS" << tag.c_str() << id; + local_histos.localPosXY[1] = td.make(); + local_histos.localPosXY[1]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Local_Position_XY_2S" << tag.c_str() << id; + local_histos.localPosXY[2] = td.make(); + local_histos.localPosXY[2]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_Mixed" << tag.c_str() << id; + local_histos.globalPosXY[0] = td.make(); + local_histos.globalPosXY[0]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_PS" << tag.c_str() << id; + local_histos.globalPosXY[1] = td.make(); + local_histos.globalPosXY[1]->SetName(histoName.str().c_str()); + + histoName.str(""); + histoName << "Global_Position_XY_2S" << tag.c_str() << id; + local_histos.globalPosXY[2] = td.make(); + local_histos.globalPosXY[2]->SetName(histoName.str().c_str()); + + /* + * Delta positions with SimHits + */ + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_Mixed" << tag.c_str() << id; + local_histos.deltaXVHSimHits[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_PS" << tag.c_str() << id; + local_histos.deltaXVHSimHits[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_2S" << tag.c_str() << id; + local_histos.deltaXVHSimHits[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_Mixed" << tag.c_str() << id; + local_histos.deltaYVHSimHits[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_PS" << tag.c_str() << id; + local_histos.deltaYVHSimHits[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_2S" << tag.c_str() << id; + local_histos.deltaYVHSimHits[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + /* + * Delta position with simHits for primary tracks only + */ + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_Mixed_P" << tag.c_str() << id; + local_histos.deltaXVHSimHits_P[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_PS_P" << tag.c_str() << id; + local_histos.deltaXVHSimHits_P[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_X_VH_SimHits_2S_P" << tag.c_str() << id; + local_histos.deltaXVHSimHits_P[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_Mixed_P" << tag.c_str() << id; + local_histos.deltaYVHSimHits_P[0] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_PS_P" << tag.c_str() << id; + local_histos.deltaYVHSimHits_P[1] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + histoName.str(""); + histoName << "Delta_Y_VH_SimHits_2S_P" << tag.c_str() << id; + local_histos.deltaYVHSimHits_P[2] = td.make(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.); + + /* + * Information on the Digis per cluster + */ + + histoName.str(""); + histoName << "Total_Digis" << tag.c_str() << id; + local_histos.totalSimHits = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + + histoName.str(""); + histoName << "Primary_Digis" << tag.c_str() << id; + local_histos.primarySimHits = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + + histoName.str(""); + histoName << "Other_Digis" << tag.c_str() << id; + local_histos.otherSimHits = td.make(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.); + + /* + * Study on the clusters combinatorial problem + */ + + histoName.str(""); + histoName << "DeltaXlocal_clusters" << tag.c_str() << id; + local_histos.deltaXlocal = td.make(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4); + histoName.str(""); + histoName << "Width" << tag.c_str() << id; + local_histos.width = td.make(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4); + histoName.str(""); + histoName << "Curvature" << tag.c_str() << id; + local_histos.curvature = td.make(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4); + + std::pair::iterator, bool> insertedIt( + histograms_.insert(std::make_pair(ival, local_histos))); + fs->file().cd("/"); + + return insertedIt.first; +} + +void VectorHitsBuilderValidation::CreateVHsXYGraph(const std::vector glVHs, + const std::vector dirVHs) { + if (glVHs.size() != dirVHs.size()) { + std::cout << "Cannot fullfil the graphs for this event. Return." << std::endl; + return; + } + + // opening canvas and drawing XY TGraph + + for (unsigned int nVH = 0; nVH < glVHs.size(); nVH++) { + //same r + if ((fabs(dirVHs.at(nVH).x()) < 10e-5) && (fabs(dirVHs.at(nVH).y()) < 10e-5)) { + continue; + + } else { + } + } + + return; +} + +void VectorHitsBuilderValidation::CreateVHsRZGraph(const std::vector glVHs, + const std::vector dirVHs) { + if (glVHs.size() != dirVHs.size()) { + std::cout << "Cannot fullfil the graphs for this event. Return." << std::endl; + return; + } + + return; +} + +void VectorHitsBuilderValidation::CreateWindowCorrGraph() { + //FIXME: This function is not working properly, yet. + + //return if we are not using Phase2 OT + if (!tkGeom_->isThere(GeomDetEnumerators::P2OTB) && !tkGeom_->isThere(GeomDetEnumerators::P2OTEC)) + return; + + for (auto det : tkGeom_->detsTOB()) { + ParallaxCorrectionRZ_->Fill(det->position().z(), det->position().perp(), 5.); + } + for (auto det : tkGeom_->detsTID()) { + ParallaxCorrectionRZ_->Fill(det->position().z(), det->position().perp(), 10.); + } + ParallaxCorrectionRZ_->Fill(0., 0., 5.); + return; +} + +unsigned int VectorHitsBuilderValidation::getLayerNumber(const DetId& detid) { + if (detid.det() == DetId::Tracker) { + if (detid.subdetId() == StripSubdetector::TOB) + return (tkTopo_->layer(detid)); + else if (detid.subdetId() == StripSubdetector::TID) + return (100 * tkTopo_->side(detid) + tkTopo_->layer(detid)); + else + return 999; + } + return 999; +} + +unsigned int VectorHitsBuilderValidation::getModuleNumber(const DetId& detid) { return (tkTopo_->module(detid)); } + +std::vector > VectorHitsBuilderValidation::getSimTrackIds( + const edm::Handle >& simLinks, const DetId& detId, uint32_t channel) const { + std::vector > simTrkId; + auto isearch = simLinks->find(detId); + if (isearch != simLinks->end()) { + // Loop over DigiSimLink in this det unit + edm::DetSet link_detset = (*isearch); + for (const auto& it : link_detset.data) { + if (channel == it.channel()) + simTrkId.push_back(std::make_pair(it.SimTrackId(), it.eventId())); + } + } + return simTrkId; +} + +unsigned int VectorHitsBuilderValidation::getSimTrackId( + const edm::Handle >& pixelSimLinks, + const DetId& detId, + unsigned int channel) const { + edm::DetSetVector::const_iterator DSViter(pixelSimLinks->find(detId)); + if (DSViter == pixelSimLinks->end()) + return 0; + for (const auto& it : DSViter->data) { + if (channel == it.channel()) + return it.SimTrackId(); + } + return 0; +} + +void VectorHitsBuilderValidation::printCluster(const GeomDetUnit* geomDetUnit, const OmniClusterRef cluster) { + if (!geomDetUnit) + return; + + const PixelGeomDetUnit* theGeomDet = dynamic_cast(geomDetUnit); + const PixelTopology& topol = theGeomDet->specificTopology(); + + unsigned int layer = getLayerNumber(geomDetUnit->geographicalId()); + unsigned int module = getModuleNumber(geomDetUnit->geographicalId()); + LogTrace("VectorHitsBuilderValidation") << "Layer:" << layer << std::endl; + if (topol.ncolumns() == 32) + LogTrace("VectorHitsBuilderValidation") << "Pixel cluster with detId:" << geomDetUnit->geographicalId().rawId() + << "(module:" << module << ") " << std::endl; + else if (topol.ncolumns() == 2) + LogTrace("VectorHitsBuilderValidation") << "Strip cluster with detId:" << geomDetUnit->geographicalId().rawId() + << "(module:" << module << ") " << std::endl; + else + std::cout << "no module?!" << std::endl; + LogTrace("VectorHitsBuilderValidation") + << "with pitch:" << topol.pitch().first << " , " << topol.pitch().second << std::endl; + LogTrace("VectorHitsBuilderValidation") << " and width:" << theGeomDet->surface().bounds().width() + << " , lenght:" << theGeomDet->surface().bounds().length() << std::endl; + + auto&& lparams = cpe_->localParameters(*cluster.cluster_phase2OT(), *theGeomDet); + + LogTrace("VectorHitsBuilderValidation") + << "\t local pos " << lparams.first << "with err " << lparams.second << std::endl; + + return; +} + +void VectorHitsBuilderValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("src", "siPhase2Clusters"); + desc.add("links", edm::InputTag("simSiPixelDigis", "Tracker")); + desc.add("VH_acc", edm::InputTag("siPhase2VectorHits", "accepted")); + desc.add("VH_rej", edm::InputTag("siPhase2VectorHits", "rejected")); + desc.add("CPE", edm::ESInputTag("phase2StripCPEESProducer", "Phase2StripCPE")); + desc.add("trackingParticleSrc", edm::InputTag("mix", "MergedTrackTruth")); + descriptions.add("vectorHitsBuilderValidation", desc); +} + +DEFINE_FWK_MODULE(VectorHitsBuilderValidation); diff --git a/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.h b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.h new file mode 100644 index 0000000000000..f0233dafe4a64 --- /dev/null +++ b/RecoLocalTracker/SiPhase2VectorHitBuilder/test/VectorHitsBuilderValidation.h @@ -0,0 +1,152 @@ +#ifndef RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitsBuilderValidation_H +#define RecoLocalTracker_SiPhase2VectorHitBuilder_VectorHitsBuilderValidation_H + +#include +#include +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h" +#include "DataFormats/SiPixelDigi/interface/PixelDigi.h" +#include "RecoLocalTracker/Phase2TrackerRecHits/interface/Phase2StripCPE.h" + +#include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" + +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "CommonTools/Utils/interface/TFileDirectory.h" + +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h" + +#include +#include +#include +#include +#include +#include +#include + +struct VHHistos { + THStack* numberVHsMixed; + TH1F* numberVHsPS; + TH1F* numberVHs2S; + + TGraph* globalPosXY[3]; + TGraph* localPosXY[3]; + + TH1F* deltaXVHSimHits[3]; + TH1F* deltaYVHSimHits[3]; + + TH1F* deltaXVHSimHits_P[3]; + TH1F* deltaYVHSimHits_P[3]; + + TH1F* digiEfficiency[3]; + + TH1F* totalSimHits; + TH1F* primarySimHits; + TH1F* otherSimHits; + + TH1F* curvature; + TH1F* width; + TH1F* deltaXlocal; +}; + +class VectorHitsBuilderValidation : public edm::one::EDAnalyzer { +public: + typedef edm::Ref, Phase2TrackerCluster1D> Phase2TrackerCluster1DRef; + + typedef std::map > SimHitsMap; + typedef std::map SimTracksMap; + + explicit VectorHitsBuilderValidation(const edm::ParameterSet&); + ~VectorHitsBuilderValidation(); + void beginJob(); + void endJob(); + void analyze(const edm::Event&, const edm::EventSetup&); + + static void fillDescriptions(edm::ConfigurationDescriptions&); + +private: + std::map::iterator createLayerHistograms(unsigned int); + void CreateVHsXYGraph(const std::vector, const std::vector); + void CreateVHsRZGraph(const std::vector, const std::vector); + void CreateWindowCorrGraph(); + + unsigned int getLayerNumber(const DetId&); + unsigned int getModuleNumber(const DetId& detid); + void printCluster(const GeomDetUnit* geomDetUnit, const OmniClusterRef cluster); + + std::pair isTrue(const VectorHit vh, + const edm::Handle >& siphase2SimLinks, + DetId& detId) const; + std::vector > getSimTrackIds( + const edm::Handle >&, const DetId&, uint32_t) const; + unsigned int getSimTrackId(const edm::Handle >& pixelSimLinks, + const DetId& detId, + unsigned int channel) const; + + edm::EDGetTokenT > srcClu_; + edm::EDGetTokenT VHacc_; + edm::EDGetTokenT VHrej_; + edm::ESInputTag cpeTag_; + const ClusterParameterEstimator* cpe_; + + edm::EDGetTokenT > siphase2OTSimLinksToken_; + edm::EDGetTokenT simHitsToken_; + edm::EDGetTokenT simTracksToken_; + edm::EDGetTokenT simVerticesToken_; + edm::EDGetTokenT trackingParticleToken_; + + const TrackerGeometry* tkGeom_; + const TrackerTopology* tkTopo_; + const MagneticField* magField_; + + TTree* tree_; + TGraph* trackerLayoutRZ_[3]; + TGraph* trackerLayoutXY_[3]; + TGraph* trackerLayoutXYBar_; + TGraph* trackerLayoutXYEC_; + TGraph* localPosXvsDeltaX_[3]; + TGraph* localPosYvsDeltaY_[3]; + TCanvas* VHXY_; + TCanvas* VHRZ_; + std::vector arrowVHs_; + + TH2D* ParallaxCorrectionRZ_; + TH1F* VHaccLayer_; + TH1F* VHrejLayer_; + TH1F* VHaccTrueLayer_; + TH1F* VHrejTrueLayer_; + TH1F* VHaccTrue_signal_Layer_; + TH1F* VHrejTrue_signal_Layer_; + + std::map histograms_; +}; +#endif diff --git a/RecoLocalTracker/SubCollectionProducers/src/SeedClusterRemoverPhase2.cc b/RecoLocalTracker/SubCollectionProducers/src/SeedClusterRemoverPhase2.cc index e8b4fe68ca477..24df700d80ed4 100644 --- a/RecoLocalTracker/SubCollectionProducers/src/SeedClusterRemoverPhase2.cc +++ b/RecoLocalTracker/SubCollectionProducers/src/SeedClusterRemoverPhase2.cc @@ -11,6 +11,7 @@ #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h" #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include "DataFormats/Common/interface/Handle.h" #include "DataFormats/Common/interface/DetSetVector.h" @@ -125,6 +126,34 @@ void SeedClusterRemoverPhase2::process(const TrackingRecHit *hit, float chi2, co assert(collectedOuterTrackers_.size() > cluster.key()); collectedOuterTrackers_[cluster.key()] = true; + } else if (hitType == typeid(VectorHit)) { + if (!doOuterTracker_) + return; + + const VectorHit *vhit = static_cast(hit); + LogDebug("SeedClusterRemoverPhase2") << "Plain VectorHit in det " << detid.rawId(); + + //lower cluster + Phase2TrackerRecHit1D::CluRef cluster = vhit->lowerCluster(); + if (cluster.id() != outerTrackerSourceProdID) + throw cms::Exception("Inconsistent Data") + << "SeedClusterRemoverPhase2: strip cluster ref from Product ID = " << cluster.id() + << " does not match with source cluster collection (ID = " << outerTrackerSourceProdID << ")\n."; + + OTs[cluster.key()] = false; + assert(collectedOuterTrackers_.size() > cluster.key()); + collectedOuterTrackers_[cluster.key()] = true; + + // upper cluster + cluster = vhit->upperCluster(); + if (cluster.id() != outerTrackerSourceProdID) + throw cms::Exception("Inconsistent Data") + << "SeedClusterRemoverPhase2: strip cluster ref from Product ID = " << cluster.id() + << " does not match with source cluster collection (ID = " << outerTrackerSourceProdID << ")\n."; + + OTs[cluster.key()] = false; + assert(collectedOuterTrackers_.size() > cluster.key()); + collectedOuterTrackers_[cluster.key()] = true; } else throw cms::Exception("NOT IMPLEMENTED") << "I received a hit that was neither SiPixelRecHit nor Phase2TrackerRecHit1D but " << hitType.name() diff --git a/RecoLocalTracker/SubCollectionProducers/src/TrackClusterRemoverPhase2.cc b/RecoLocalTracker/SubCollectionProducers/src/TrackClusterRemoverPhase2.cc index a7afaf340eba6..18e356a0520b8 100644 --- a/RecoLocalTracker/SubCollectionProducers/src/TrackClusterRemoverPhase2.cc +++ b/RecoLocalTracker/SubCollectionProducers/src/TrackClusterRemoverPhase2.cc @@ -18,6 +18,7 @@ #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackerRecHit2D/interface/ClusterRemovalInfo.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include "TrackingTools/PatternTools/interface/TrackCollectionTokens.h" @@ -190,6 +191,22 @@ namespace { collectedPixels[cluster.key()] = true; else if (cluster.isPhase2()) collectedPhase2OTs[cluster.key()] = true; + + // Phase 2 OT is defined as Pixel detector (for now) + const auto& hitType = typeid(hit); + if (hitType == typeid(VectorHit)) { + auto const& vectorHit = reinterpret_cast(hit); + auto const& lowCluster = vectorHit.lowerClusterRef(); + auto const& uppCluster = vectorHit.upperClusterRef(); + LogTrace("TrackClusterRemoverPhase2") + << "masking a VHit with lowCluster key: " << lowCluster.key() << " and upper key: " << uppCluster.key(); + if (lowCluster.isPhase2()) + collectedPhase2OTs[lowCluster.key()] = true; + if (uppCluster.isPhase2()) + collectedPhase2OTs[uppCluster.key()] = true; + } else { + LogTrace("TrackClusterRemoverPhase2") << "it is not a VHit."; + } } } diff --git a/RecoTracker/CkfPattern/src/PrintoutHelper.cc b/RecoTracker/CkfPattern/src/PrintoutHelper.cc index 807f90e009355..bae0056c13a86 100644 --- a/RecoTracker/CkfPattern/src/PrintoutHelper.cc +++ b/RecoTracker/CkfPattern/src/PrintoutHelper.cc @@ -56,9 +56,7 @@ std::string PrintoutHelper::dumpMeasurement(const TrajectoryMeasurement& tm) { << "p: " << tm.predictedState().globalMomentum() << "\n"; else buffer << "no valid state\n"; - buffer - // <<"geomdet pointer from rechit: "<det()<<"\n" - << "detId: " << tm.recHit()->geographicalId().rawId(); + buffer << "detId: " << tm.recHit()->geographicalId().rawId(); if (tm.recHit()->isValid()) { buffer << "\n hit global x: " << tm.recHit()->globalPosition() << "\n hit global error: " << tm.recHit()->globalPositionError().matrix() @@ -75,19 +73,6 @@ std::string PrintoutHelper::dumpMeasurement(const TrajectoryMeasurement& tm) { std::string PrintoutHelper::regressionTest(const TrackerGeometry& tracker, std::vector& unsmoothedResult) { std::stringstream buffer; - /* - for(iseed=theSeedColl.begin();iseed!=theSeedColl.end();iseed++){ - DetId tmpId = DetId( iseed->startingState().detId()); - const GeomDet* tmpDet = tracker->idToDet( tmpId ); - GlobalVector gv = tmpDet->surface().toGlobal( iseed->startingState().parameters().momentum() ); - - LogTrace("TrackingRegressionTest") << "seed perp,phi,eta : " - << gv.perp() << " , " - << gv.phi() << " , " - << gv.eta() ; - } - */ - buffer << "number of finalTrajectories: " << unsmoothedResult.size() << std::endl; for (std::vector::const_iterator it = unsmoothedResult.begin(); it != unsmoothedResult.end(); it++) { if (it->lastMeasurement().updatedState().isValid()) { @@ -102,5 +87,12 @@ std::string PrintoutHelper::regressionTest(const TrackerGeometry& tracker, std:: buffer << "candidate with invalid last measurement state!" << std::endl; } buffer << "================================================="; + buffer << "=========== Traj in details =====================\n"; + for (const auto& it : unsmoothedResult) { + for (const auto& hit : it.measurements()) { + buffer << "measurement : " << hit.recHit()->geographicalId().rawId() << std::endl; + } + buffer << "================\n"; + } return buffer.str(); } diff --git a/RecoTracker/FinalTrackSelectors/python/earlyGeneralTracks_cfi.py b/RecoTracker/FinalTrackSelectors/python/earlyGeneralTracks_cfi.py index 6ab340983ab9c..294008807c0df 100644 --- a/RecoTracker/FinalTrackSelectors/python/earlyGeneralTracks_cfi.py +++ b/RecoTracker/FinalTrackSelectors/python/earlyGeneralTracks_cfi.py @@ -101,3 +101,12 @@ makeReKeyedSeeds = cms.untracked.bool(False) ) ) +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +def _extend_pixelLess(x): + x.TrackProducers += ['pixelLessStepTracks'] + x.hasSelector += [1] + x.indivShareFrac += [0.095] + x.selectedTrackQuals += ['pixelLessStepSelector:pixelLessStep'] + x.setsToMerge[0].tLists += [6] +(trackingPhase2PU140 & vectorHits).toModify(earlyGeneralTracks, _extend_pixelLess) + diff --git a/RecoTracker/GeometryESProducer/plugins/TrackerRecoGeometryESProducer.cc b/RecoTracker/GeometryESProducer/plugins/TrackerRecoGeometryESProducer.cc index ee6462380a982..04f9f7fba6026 100644 --- a/RecoTracker/GeometryESProducer/plugins/TrackerRecoGeometryESProducer.cc +++ b/RecoTracker/GeometryESProducer/plugins/TrackerRecoGeometryESProducer.cc @@ -28,11 +28,13 @@ class TrackerRecoGeometryESProducer : public edm::ESProducer { private: edm::ESGetToken geomToken_; edm::ESGetToken tTopToken_; + bool usePhase2Stacks_; }; using namespace edm; -TrackerRecoGeometryESProducer::TrackerRecoGeometryESProducer(const edm::ParameterSet &p) { +TrackerRecoGeometryESProducer::TrackerRecoGeometryESProducer(const edm::ParameterSet &p) + : usePhase2Stacks_(p.getParameter("usePhase2Stacks")) { auto c = setWhatProduced(this); // 08-Oct-2007 - Patrick Janot @@ -50,12 +52,14 @@ std::unique_ptr TrackerRecoGeometryESProducer::produce( TrackerGeometry const &tG = iRecord.get(geomToken_); GeometricSearchTrackerBuilder builder; - return std::unique_ptr(builder.build(tG.trackerDet(), &tG, &iRecord.get(tTopToken_))); + return std::unique_ptr( + builder.build(tG.trackerDet(), &tG, &iRecord.get(tTopToken_), usePhase2Stacks_)); } void TrackerRecoGeometryESProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { edm::ParameterSetDescription desc; + desc.add("usePhase2Stacks", false); desc.addUntracked("trackerGeometryLabel", ""); descriptions.addDefault(desc); } diff --git a/RecoTracker/GeometryESProducer/python/TrackerRecoGeometryESProducer_cfi.py b/RecoTracker/GeometryESProducer/python/TrackerRecoGeometryESProducer_cfi.py index fb33c1fff2fc2..048fd49b6da54 100644 --- a/RecoTracker/GeometryESProducer/python/TrackerRecoGeometryESProducer_cfi.py +++ b/RecoTracker/GeometryESProducer/python/TrackerRecoGeometryESProducer_cfi.py @@ -1,5 +1,9 @@ import FWCore.ParameterSet.Config as cms -TrackerRecoGeometryESProducer = cms.ESProducer("TrackerRecoGeometryESProducer") +TrackerRecoGeometryESProducer = cms.ESProducer("TrackerRecoGeometryESProducer", + usePhase2Stacks = cms.bool(False) +) +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(TrackerRecoGeometryESProducer, usePhase2Stacks = True) diff --git a/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py b/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py index dad46627f08d7..ad72a25b1959c 100644 --- a/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py +++ b/RecoTracker/IterativeTracking/python/HighPtTripletStep_cff.py @@ -318,6 +318,9 @@ ] #end of vpset ) #end of clone +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(highPtTripletStepSelector.trackSelectors[2], minNumberLayers = 3, minNumber3DLayers = 3, d0_par1 = ( 0.5, 4.0 ), dz_par1 = ( 0.6, 4.0 )) + # Final sequence HighPtTripletStepTask = cms.Task(highPtTripletStepClusters, highPtTripletStepSeedLayers, @@ -327,7 +330,6 @@ highPtTripletStepSeeds, highPtTripletStepTrackCandidates, highPtTripletStepTracks, -# highPtTripletStepClassifier1,highPtTripletStepClassifier2,highPtTripletStepClassifier3* highPtTripletStep) HighPtTripletStep = cms.Sequence(HighPtTripletStepTask) _HighPtTripletStepTask_Phase2PU140 = HighPtTripletStepTask.copy() diff --git a/RecoTracker/IterativeTracking/python/LowPtTripletStep_cff.py b/RecoTracker/IterativeTracking/python/LowPtTripletStep_cff.py index fdd1f29f75afe..266a80b156123 100644 --- a/RecoTracker/IterativeTracking/python/LowPtTripletStep_cff.py +++ b/RecoTracker/IterativeTracking/python/LowPtTripletStep_cff.py @@ -363,6 +363,9 @@ ) #end of clone +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(lowPtTripletStepSelector.trackSelectors[2], minNumberLayers = 3, minNumber3DLayers = 3) + # Final sequence LowPtTripletStepTask = cms.Task(lowPtTripletStepClusters, diff --git a/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py b/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py index 4180e142dd8ff..904b1d81aedf0 100644 --- a/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py +++ b/RecoTracker/IterativeTracking/python/PixelLessStep_cff.py @@ -104,7 +104,30 @@ MTID = None, MTEC = None, ) - +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(pixelLessStepSeedLayers, + layerList = [ + 'TOB1+TOB2', 'TOB2+TOB3', +# 'TOB3+TOB4', 'TOB4+TOB5', + 'TID1_pos+TID2_pos', 'TID1_neg+TID2_neg' + ], + TOB = cms.PSet( + TTRHBuilder = cms.string('WithTrackAngle'), + clusterChargeCut = cms.PSet(refToPSet_ = cms.string('SiStripClusterChargeCutNone')), + vectorRecHits = cms.InputTag("siPhase2VectorHits", 'vectorHitsAccepted'), + skipClusters = cms.InputTag('pixelLessStepClusters') + ), + TIB = None, + TID = dict( + clusterChargeCut = dict(refToPSet_ = cms.string('SiStripClusterChargeCutNone')), + vectorRecHits = cms.InputTag("siPhase2VectorHits", 'accepted'), + maxRing = cms.int32(8) + ), + TEC = None, + MTIB = None, + MTID = None, + MTEC = None, +) # TrackingRegion from RecoTracker.TkTrackingRegions.globalTrackingRegionFromBeamSpotFixedZ_cfi import globalTrackingRegionFromBeamSpotFixedZ as _globalTrackingRegionFromBeamSpotFixedZ pixelLessStepTrackingRegions = _globalTrackingRegionFromBeamSpotFixedZ.clone(RegionPSet = dict( @@ -194,6 +217,16 @@ )) fastSim.toReplaceWith(pixelLessStepSeeds,_fastSim_pixelLessStepSeeds) +vectorHits.toModify(pixelLessStepHitDoublets, produceSeedingHitSets=True, produceIntermediateHitDoublets=False) +vectorHits.toModify(pixelLessStepSeeds, + seedingHitSets = "pixelLessStepHitDoublets", + SeedComparitorPSet = dict( + ClusterShapeHitFilterName = cms.string('ClusterShapeHitFilter'), + FilterAtHelixStage = cms.bool(False), + FilterStripHits = cms.bool(False), + ) +) + # QUALITY CUTS DURING TRACK BUILDING import TrackingTools.TrajectoryFiltering.TrajectoryFilter_cff _pixelLessStepTrajectoryFilterBase = TrackingTools.TrajectoryFiltering.TrajectoryFilter_cff.CkfBaseTrajectoryFilter_block.clone( @@ -208,6 +241,9 @@ for e in [pp_on_XeXe_2017, pp_on_AA_2018]: e.toModify(pixelLessStepTrajectoryFilter, minPt=2.0) +vectorHits.toReplaceWith(pixelLessStepTrajectoryFilter, _pixelLessStepTrajectoryFilterBase) +vectorHits.toModify(pixelLessStepTrajectoryFilter, minimumNumberOfHits = 4, maxLostHits = 1) + import RecoTracker.MeasurementDet.Chi2ChargeMeasurementEstimator_cfi pixelLessStepChi2Est = RecoTracker.MeasurementDet.Chi2ChargeMeasurementEstimator_cfi.Chi2ChargeMeasurementEstimator.clone( ComponentName = 'pixelLessStepChi2Est', @@ -219,6 +255,11 @@ clusterChargeCut = dict(refToPSet_ = 'SiStripClusterChargeCutTiny') ) +vectorHits.toModify(pixelLessStepChi2Est, + clusterChargeCut = dict(refToPSet_ = 'SiStripClusterChargeCutNone'), + MaxChi2 = 30.0 +) + # TRACK BUILDING import RecoTracker.CkfPattern.GroupedCkfTrajectoryBuilder_cfi pixelLessStepTrajectoryBuilder = RecoTracker.CkfPattern.GroupedCkfTrajectoryBuilder_cfi.GroupedCkfTrajectoryBuilder.clone( @@ -252,6 +293,11 @@ ) ) +vectorHits.toModify(pixelLessStepTrackCandidates, + phase2clustersToSkip = cms.InputTag('pixelLessStepClusters'), + clustersToSkip = None +) + from TrackingTools.TrajectoryCleaning.TrajectoryCleanerBySharedHits_cfi import trajectoryCleanerBySharedHits pixelLessStepTrajectoryCleanerBySharedHits = trajectoryCleanerBySharedHits.clone( ComponentName = 'pixelLessStepTrajectoryCleanerBySharedHits', @@ -358,6 +404,55 @@ vertices = 'pixelVertices'#end of vpset ) #end of clone +vectorHits.toModify(pixelLessStepSelector, + GBRForestLabel = None, + useAnyMVA = None, + trackSelectors= cms.VPSet( + RecoTracker.FinalTrackSelectors.multiTrackSelector_cfi.looseMTS.clone( + name = 'pixelLessStepLoose', + chi2n_par = 1.0, + res_par = ( 0.003, 0.001 ), + minNumberLayers = 0, + maxNumberLostLayers = 1, + minNumber3DLayers = 0, + d0_par1 = ( 0.9, 4.0 ), + dz_par1 = ( 0.9, 4.0 ), + d0_par2 = ( 1.0, 4.0 ), + dz_par2 = ( 1.0, 4.0 ) + ), + RecoTracker.FinalTrackSelectors.multiTrackSelector_cfi.tightMTS.clone( + name = 'pixelLessStepTight', + preFilterName = 'pixelLessStepLoose', + chi2n_par = 0.35, + res_par = ( 0.003, 0.001 ), + minNumberLayers = 4, + maxNumberLostLayers = 0, + minNumber3DLayers = 3, + d0_par1 = ( 1.1, 4.0 ), + dz_par1 = ( 1.1, 4.0 ), + d0_par2 = ( 1.1, 4.0 ), + dz_par2 = ( 1.1, 4.0 ) + ), + RecoTracker.FinalTrackSelectors.multiTrackSelector_cfi.highpurityMTS.clone( + name = 'QualityMasks', + preFilterName = 'pixelLessStepTight', + chi2n_par = 0.2, + res_par = ( 0.003, 0.001 ), + minNumberLayers = 1, + maxNumberLostLayers = 2, + minNumber3DLayers = 0, + d0_par1 = ( 100., 4.0 ), + dz_par1 = ( 100., 4.0 ), + d0_par2 = ( 100., 4.0 ), + dz_par2 = ( 100., 4.0 ) + ), + ), + vertices = 'firstStepPrimaryVertices' +) + +vectorHits.toModify(pixelLessStepSelector.trackSelectors[2], name = 'pixelLessStep') + + PixelLessStepTask = cms.Task(pixelLessStepClusters, pixelLessStepSeedLayers, pixelLessStepTrackingRegions, @@ -386,3 +481,17 @@ ,pixelLessStep ) ) +from RecoLocalTracker.SiPhase2VectorHitBuilder.siPhase2VectorHits_cfi import * +from RecoTracker.TkSeedGenerator.SeedingOTEDProducer_cfi import SeedingOTEDProducer as _SeedingOTEDProducer +pixelLessStepSeeds_vectorHits = _SeedingOTEDProducer.clone() + +_PixelLessStepTask_vectorHits = cms.Task(siPhase2VectorHits, + pixelLessStepClusters, + pixelLessStepSeeds, + pixelLessStepTrackCandidates, + pixelLessStepTracks, + pixelLessStepSelector) +_PixelLessStep_vectorHits = cms.Sequence(_PixelLessStepTask_vectorHits) +vectorHits.toReplaceWith(pixelLessStepSeeds,pixelLessStepSeeds_vectorHits) +vectorHits.toReplaceWith(PixelLessStepTask, _PixelLessStepTask_vectorHits) +vectorHits.toReplaceWith(PixelLessStep, _PixelLessStep_vectorHits) diff --git a/RecoTracker/IterativeTracking/python/PixelPairStep_cff.py b/RecoTracker/IterativeTracking/python/PixelPairStep_cff.py index bdca6fd4af1e7..81eb716d6a39b 100644 --- a/RecoTracker/IterativeTracking/python/PixelPairStep_cff.py +++ b/RecoTracker/IterativeTracking/python/PixelPairStep_cff.py @@ -419,6 +419,8 @@ vertices = 'firstStepPrimaryVertices' ) #end of clone +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(pixelPairStepSelector.trackSelectors[2], minNumberLayers = 3, minNumber3DLayers = 3) # Final sequence PixelPairStepTask = cms.Task(pixelPairStepClusters, diff --git a/RecoTracker/IterativeTracking/python/iterativeTkConfig.py b/RecoTracker/IterativeTracking/python/iterativeTkConfig.py index a34b00031488e..0f187dcf4184b 100644 --- a/RecoTracker/IterativeTracking/python/iterativeTkConfig.py +++ b/RecoTracker/IterativeTracking/python/iterativeTkConfig.py @@ -56,6 +56,8 @@ "DetachedQuadStep", "PixelPairStep", ] +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(_iterations_trackingPhase2PU140, func=lambda x: x.append('PixelLessStep')) _iterations_muonSeeded = [ "MuonSeededStepInOut", "MuonSeededStepOutIn", diff --git a/RecoTracker/MeasurementDet/BuildFile.xml b/RecoTracker/MeasurementDet/BuildFile.xml index 5ed4491d1708d..a7ea54e101d9b 100644 --- a/RecoTracker/MeasurementDet/BuildFile.xml +++ b/RecoTracker/MeasurementDet/BuildFile.xml @@ -22,3 +22,4 @@ + diff --git a/RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h b/RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h index c65203dba5624..a7cf35cbab0cb 100644 --- a/RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h +++ b/RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h @@ -8,6 +8,7 @@ class Phase2OTMeasurementDetSet; #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include "DataFormats/Common/interface/ContainerMask.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" @@ -25,6 +26,8 @@ class MeasurementTrackerEvent { const StMeasurementDetSet *strips, const PxMeasurementDetSet *pixels, const Phase2OTMeasurementDetSet *phase2OT, + const VectorHitCollection *phase2OTVectorHits, + const VectorHitCollection *phase2OTVectorHitsRej, const std::vector &stripClustersToSkip, const std::vector &pixelClustersToSkip, const std::vector &phase2OTClustersToSkip) @@ -32,6 +35,8 @@ class MeasurementTrackerEvent { theStripData(strips), thePixelData(pixels), thePhase2OTData(phase2OT), + thePhase2OTVectorHits(phase2OTVectorHits), + thePhase2OTVectorHitsRej(phase2OTVectorHitsRej), theOwner(true), theStripClustersToSkip(stripClustersToSkip), thePixelClustersToSkip(pixelClustersToSkip), @@ -57,6 +62,8 @@ class MeasurementTrackerEvent { const StMeasurementDetSet &stripData() const { return *theStripData; } const PxMeasurementDetSet &pixelData() const { return *thePixelData; } const Phase2OTMeasurementDetSet &phase2OTData() const { return *thePhase2OTData; } + const VectorHitCollection &phase2OTVectorHits() const { return *thePhase2OTVectorHits; } + const VectorHitCollection &phase2OTVectorHitsRej() const { return *thePhase2OTVectorHitsRej; } const std::vector &stripClustersToSkip() const { return theStripClustersToSkip; } const std::vector &pixelClustersToSkip() const { return thePixelClustersToSkip; } const std::vector &phase2OTClustersToSkip() const { return thePhase2OTClustersToSkip; } @@ -73,6 +80,8 @@ class MeasurementTrackerEvent { const StMeasurementDetSet *theStripData = nullptr; const PxMeasurementDetSet *thePixelData = nullptr; const Phase2OTMeasurementDetSet *thePhase2OTData = nullptr; + const VectorHitCollection *thePhase2OTVectorHits = nullptr; + const VectorHitCollection *thePhase2OTVectorHitsRej = nullptr; bool theOwner = false; // do I own the tree above? // these could be const pointers as well, but ContainerMask doesn't expose the vector std::vector theStripClustersToSkip; diff --git a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerESProducer.cc b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerESProducer.cc index 9087a115fa534..54b75367b7df5 100644 --- a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerESProducer.cc +++ b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerESProducer.cc @@ -10,6 +10,7 @@ #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" #include "RecoLocalTracker/Phase2TrackerRecHits/interface/Phase2StripCPE.h" #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripRecHitMatcher.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" #include "RecoTracker/Record/interface/CkfComponentsRecord.h" diff --git a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.cc b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.cc index d8a3b89ae5c97..29b49c425832d 100644 --- a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.cc +++ b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.cc @@ -39,7 +39,8 @@ MeasurementTrackerEventProducer::MeasurementTrackerEventProducer(const edm::Para edm::InputTag skip = iConfig.getParameter("skipClusters"); selfUpdateSkipClusters_ = !(skip == edm::InputTag("")); LogDebug("MeasurementTracker") << "skipping clusters: " << selfUpdateSkipClusters_; - isPhase2 = false; + isPhase2_ = false; + useVectorHits_ = false; if (!iConfig.getParameter("stripClusterProducer").empty()) { theStripClusterLabel = consumes>( @@ -58,7 +59,14 @@ MeasurementTrackerEventProducer::MeasurementTrackerEventProducer(const edm::Para if (!iConfig.getParameter("Phase2TrackerCluster1DProducer").empty()) { thePh2OTClusterLabel = consumes>( edm::InputTag(iConfig.getParameter("Phase2TrackerCluster1DProducer"))); - isPhase2 = true; + isPhase2_ = true; + } + if (!(iConfig.getParameter("vectorHits") == edm::InputTag("") || + iConfig.getParameter("vectorHitsRej") == edm::InputTag(""))) { + thePh2OTVectorHitsLabel = consumes(iConfig.getParameter("vectorHits")); + thePh2OTVectorHitsRejLabel = consumes(iConfig.getParameter("vectorHitsRej")); + isPhase2_ = true; + useVectorHits_ = true; } produces(); @@ -72,6 +80,8 @@ void MeasurementTrackerEventProducer::fillDescriptions(edm::ConfigurationDescrip desc.add("pixelClusterProducer", "siPixelClusters"); desc.add("stripClusterProducer", "siStripClusters"); desc.add("Phase2TrackerCluster1DProducer", ""); + desc.add("vectorHits", edm::InputTag("")); + desc.add("vectorHitsRej", edm::InputTag("")); desc.add>("inactivePixelDetectorLabels", std::vector{{edm::InputTag("siPixelDigis")}}) @@ -97,6 +107,7 @@ void MeasurementTrackerEventProducer::produce(edm::Event& iEvent, const edm::Eve auto stripData = std::make_unique(measurementTracker->stripDetConditions()); auto pixelData = std::make_unique(measurementTracker->pixelDetConditions()); auto phase2OTData = std::make_unique(measurementTracker->phase2DetConditions()); + std::vector stripClustersToSkip; std::vector pixelClustersToSkip; std::vector phase2ClustersToSkip; @@ -112,10 +123,16 @@ void MeasurementTrackerEventProducer::produce(edm::Event& iEvent, const edm::Eve // put into MTE // put into event + // + + const VectorHitCollection* phase2OTVectorHits = useVectorHits_ ? &iEvent.get(thePh2OTVectorHitsLabel) : nullptr; + const VectorHitCollection* phase2OTVectorHitsRej = useVectorHits_ ? &iEvent.get(thePh2OTVectorHitsRejLabel) : nullptr; iEvent.put(std::make_unique(*measurementTracker, stripData.release(), pixelData.release(), phase2OTData.release(), + phase2OTVectorHits, + phase2OTVectorHitsRej, stripClustersToSkip, pixelClustersToSkip, phase2ClustersToSkip)); @@ -350,8 +367,10 @@ void MeasurementTrackerEventProducer::updateStrips(const edm::Event& event, //FIXME: just a temporary solution for phase2! void MeasurementTrackerEventProducer::updatePhase2OT(const edm::Event& event, Phase2OTMeasurementDetSet& thePh2OTDets) const { + thePh2OTDets.setEmpty(); + // Phase2OT Clusters - if (isPhase2) { + if (isPhase2_) { if (thePh2OTClusterLabel.isUninitialized()) { //clusters have not been produced thePh2OTDets.setActiveThisEvent(false); } else { diff --git a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.h b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.h index c502477a8e3d8..53ee67905c456 100644 --- a/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.h +++ b/RecoTracker/MeasurementDet/plugins/MeasurementTrackerEventProducer.h @@ -12,6 +12,7 @@ #include "DataFormats/DetId/interface/DetIdCollection.h" #include "DataFormats/SiPixelDetId/interface/PixelFEDChannel.h" #include "CondFormats/SiPixelObjects/interface/SiPixelFedCablingMap.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" class dso_hidden MeasurementTrackerEventProducer final : public edm::stream::EDProducer<> { public: @@ -40,6 +41,8 @@ class dso_hidden MeasurementTrackerEventProducer final : public edm::stream::EDP edm::EDGetTokenT> thePixelClusterLabel; edm::EDGetTokenT> theStripClusterLabel; edm::EDGetTokenT> thePh2OTClusterLabel; + edm::EDGetTokenT thePh2OTVectorHitsLabel; + edm::EDGetTokenT thePh2OTVectorHitsRejLabel; edm::EDGetTokenT>> thePixelClusterMask; edm::EDGetTokenT>> theStripClusterMask; @@ -50,7 +53,8 @@ class dso_hidden MeasurementTrackerEventProducer final : public edm::stream::EDP bool selfUpdateSkipClusters_; bool switchOffPixelsIfEmpty_; - bool isPhase2; + bool isPhase2_; + bool useVectorHits_; }; #endif diff --git a/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.cc b/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.cc index aab7c42efca4c..f284d0c56d30f 100644 --- a/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.cc +++ b/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.cc @@ -1,17 +1,18 @@ #include "TkStackMeasurementDet.h" #include "TrackingTools/MeasurementDet/interface/MeasurementDetException.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" using namespace std; TkStackMeasurementDet::TkStackMeasurementDet(const StackGeomDet* gdet, const PixelClusterParameterEstimator* cpe) - : MeasurementDet(gdet), thePixelCPE(cpe), theInnerDet(nullptr), theOuterDet(nullptr) {} + : MeasurementDet(gdet), thePixelCPE(cpe), theLowerDet(nullptr), theUpperDet(nullptr) {} void TkStackMeasurementDet::init(const MeasurementDet* lowerDet, const MeasurementDet* upperDet) { - theInnerDet = dynamic_cast(lowerDet); - theOuterDet = dynamic_cast(upperDet); + theLowerDet = dynamic_cast(lowerDet); + theUpperDet = dynamic_cast(upperDet); - if ((theInnerDet == nullptr) || (theOuterDet == nullptr)) { + if ((theLowerDet == nullptr) || (theUpperDet == nullptr)) { throw MeasurementDetException( "TkStackMeasurementDet ERROR: Trying to glue a det which is not a TkPhase2OTMeasurementDet"); } @@ -20,10 +21,80 @@ void TkStackMeasurementDet::init(const MeasurementDet* lowerDet, const Measureme TkStackMeasurementDet::RecHitContainer TkStackMeasurementDet::recHits(const TrajectoryStateOnSurface& ts, const MeasurementTrackerEvent& data) const { RecHitContainer result; - /* - HitCollectorForRecHits collector( &fastGeomDet(), theMatcher, theCPE, result ); - collectRecHits(ts, collector); -*/ + + if (data.phase2OTVectorHits().empty()) + return result; + LogTrace("MeasurementTracker") << " is not empty"; + if (!isActive(data)) + return result; + LogTrace("MeasurementTracker") << " and is active"; + + //find clusters to skip + const detset& lowerDetSet = data.phase2OTData().detSet(lowerDet()->index()); + const detset& upperDetSet = data.phase2OTData().detSet(upperDet()->index()); + std::vector skipClustersUpper(data.phase2OTClustersToSkip().empty() ? 0 : upperDetSet.size(), false); + std::vector skipClustersLower(data.phase2OTClustersToSkip().empty() ? 0 : lowerDetSet.size(), false); + + const Phase2TrackerCluster1D* begin = nullptr; + if (!data.phase2OTData().handle()->data().empty()) { + begin = &(data.phase2OTData().handle()->data().front()); + } + if (!data.phase2OTClustersToSkip().empty()) { + if (!lowerDetSet.empty()) { + for (const_iterator cil = lowerDetSet.begin(); cil != lowerDetSet.end(); ++cil) { + if (cil < begin) { + edm::LogError("IndexMisMatch") << "TkStackMeasurementDet cannot create hit because of index mismatch."; + return result; + } + unsigned int indexl = cil - begin; + if (data.phase2OTClustersToSkip()[indexl]) { + int iLocalL = std::distance(lowerDetSet.begin(), cil); + skipClustersLower[iLocalL] = true; + } + } + } + if (!upperDetSet.empty()) { + for (const_iterator ciu = upperDetSet.begin(); ciu != upperDetSet.end(); ++ciu) { + if (ciu < begin) { + edm::LogError("IndexMisMatch") << "TkStackMeasurementDet cannot create hit because of index mismatch."; + return result; + } + unsigned int indexu = ciu - begin; + if (data.phase2OTClustersToSkip()[indexu]) { + int iLocalU = std::distance(upperDetSet.begin(), ciu); + skipClustersUpper[iLocalU] = true; + } + } + } + } + DetId detIdStack = specificGeomDet().geographicalId(); + + auto iterator = data.phase2OTVectorHits().find(detIdStack); + if (iterator == data.phase2OTVectorHits().end()) + return result; + for (const auto& vecHit : data.phase2OTVectorHits()[detIdStack]) { + if (!data.phase2OTClustersToSkip().empty()) { + if (skipClustersLower[vecHit.lowerCluster().key() - lowerDetSet.offset()]) + continue; + if (skipClustersUpper[vecHit.upperCluster().key() - upperDetSet.offset()]) + continue; + } + result.push_back(std::make_shared(vecHit)); + } + + iterator = data.phase2OTVectorHitsRej().find(detIdStack); + if (iterator == data.phase2OTVectorHitsRej().end()) + return result; + for (const auto& vecHit : data.phase2OTVectorHitsRej()[detIdStack]) { + if (!data.phase2OTClustersToSkip().empty()) { + if (skipClustersLower[vecHit.lowerCluster().key() - lowerDetSet.offset()]) + continue; + if (skipClustersUpper[vecHit.upperCluster().key() - upperDetSet.offset()]) + continue; + } + result.push_back(std::make_shared(vecHit)); + } + return result; } @@ -31,5 +102,31 @@ bool TkStackMeasurementDet::measurements(const TrajectoryStateOnSurface& stateOn const MeasurementEstimator& est, const MeasurementTrackerEvent& data, TempMeasurements& result) const { - return true; + LogDebug("MeasurementTracker") << "TkStackMeasurementDet::measurements"; + + if (!isActive(data)) { + result.add(theInactiveHit, 0.F); + return true; + } + + LogTrace("MeasurementTracker") << " is active"; + + auto oldSize = result.size(); + MeasurementDet::RecHitContainer&& allHits = recHits(stateOnThisDet, data); + + for (auto&& hit : allHits) { + std::pair diffEst = est.estimate(stateOnThisDet, *hit); + if (diffEst.first) { + LogDebug("MeasurementTracker") << "New vh added with chi2: " << diffEst.second; + result.add(std::move(hit), diffEst.second); + } + } + + if (result.size() > oldSize) + return true; + + // create a TrajectoryMeasurement with an invalid RecHit and zero estimate + result.add(theMissingHit, 0.F); + LogDebug("MeasurementTracker") << "adding missing hit"; + return false; } diff --git a/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.h b/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.h index 91129ab61732c..7dc973344b529 100644 --- a/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.h +++ b/RecoTracker/MeasurementDet/plugins/TkStackMeasurementDet.h @@ -5,6 +5,7 @@ #include "TkPhase2OTMeasurementDet.h" #include "Geometry/CommonDetUnit/interface/StackGeomDet.h" +#include "RecoLocalTracker/SiPhase2VectorHitBuilder/interface/VectorHitBuilderAlgorithm.h" #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" #include "FWCore/Utilities/interface/Visibility.h" @@ -19,19 +20,26 @@ class TkStackMeasurementDet final : public MeasurementDet { RecHitContainer recHits(const TrajectoryStateOnSurface&, const MeasurementTrackerEvent& data) const override; const StackGeomDet& specificGeomDet() const { return static_cast(fastGeomDet()); } + typedef edm::Ref, Phase2TrackerCluster1D> Phase2TrackerCluster1DRef; + + typedef edmNew::DetSet detset; + typedef detset::const_iterator const_iterator; bool measurements(const TrajectoryStateOnSurface& stateOnThisDet, const MeasurementEstimator& est, const MeasurementTrackerEvent& data, TempMeasurements& result) const override; - const TkPhase2OTMeasurementDet* lowerDet() const { return theInnerDet; } - const TkPhase2OTMeasurementDet* upperDet() const { return theOuterDet; } + const TkPhase2OTMeasurementDet* lowerDet() const { return theLowerDet; } + const TkPhase2OTMeasurementDet* upperDet() const { return theUpperDet; } /// return TRUE if both lower and upper components are active bool isActive(const MeasurementTrackerEvent& data) const override { return lowerDet()->isActive(data) && upperDet()->isActive(data); } + bool isEmpty(const Phase2OTMeasurementDetSet& data) const { + return data.empty(lowerDet()->index()) || data.empty(upperDet()->index()); + } /// return TRUE if at least one of the lower and upper components has badChannels bool hasBadComponents(const TrajectoryStateOnSurface& tsos, const MeasurementTrackerEvent& data) const override { @@ -40,8 +48,8 @@ class TkStackMeasurementDet final : public MeasurementDet { private: const PixelClusterParameterEstimator* thePixelCPE; - const TkPhase2OTMeasurementDet* theInnerDet; - const TkPhase2OTMeasurementDet* theOuterDet; + const TkPhase2OTMeasurementDet* theLowerDet; + const TkPhase2OTMeasurementDet* theUpperDet; }; #endif diff --git a/RecoTracker/MeasurementDet/python/MeasurementTrackerEventProducer_cfi.py b/RecoTracker/MeasurementDet/python/MeasurementTrackerEventProducer_cfi.py index 6442eb1de605a..2d4e807b21872 100644 --- a/RecoTracker/MeasurementDet/python/MeasurementTrackerEventProducer_cfi.py +++ b/RecoTracker/MeasurementDet/python/MeasurementTrackerEventProducer_cfi.py @@ -23,7 +23,11 @@ inactiveStripDetectorLabels = cms.VInputTag(), switchOffPixelsIfEmpty = False ) - +from Configuration.ProcessModifiers.vectorHits_cff import vectorHits +vectorHits.toModify(MeasurementTrackerEvent, + vectorHits = ("siPhase2VectorHits", "accepted"), + vectorHitsRej = ("siPhase2VectorHits", "rejected"), +) MeasurementTrackerEventPreSplitting = MeasurementTrackerEvent.clone( pixelClusterProducer = 'siPixelClustersPreSplitting' diff --git a/RecoTracker/MeasurementDet/src/MeasurementTrackerEvent.cc b/RecoTracker/MeasurementDet/src/MeasurementTrackerEvent.cc index 274a2f62cfbde..e764ff3607fc7 100644 --- a/RecoTracker/MeasurementDet/src/MeasurementTrackerEvent.cc +++ b/RecoTracker/MeasurementDet/src/MeasurementTrackerEvent.cc @@ -18,6 +18,8 @@ MeasurementTrackerEvent::MeasurementTrackerEvent(MeasurementTrackerEvent &&other theStripData = other.theStripData; thePixelData = other.thePixelData; thePhase2OTData = other.thePhase2OTData; + thePhase2OTVectorHits = other.thePhase2OTVectorHits; + thePhase2OTVectorHitsRej = other.thePhase2OTVectorHitsRej; theOwner = other.theOwner; other.theOwner = false; // make sure to fully transfer the ownership theStripClustersToSkip = std::move(other.theStripClustersToSkip); @@ -28,6 +30,8 @@ MeasurementTrackerEvent &MeasurementTrackerEvent::operator=(MeasurementTrackerEv theStripData = other.theStripData; thePixelData = other.thePixelData; thePhase2OTData = other.thePhase2OTData; + thePhase2OTVectorHits = other.thePhase2OTVectorHits; + thePhase2OTVectorHitsRej = other.thePhase2OTVectorHitsRej; theOwner = other.theOwner; other.theOwner = false; // make sure to fully transfer the ownership theStripClustersToSkip = std::move(other.theStripClustersToSkip); @@ -43,6 +47,8 @@ MeasurementTrackerEvent::MeasurementTrackerEvent( theStripData(trackerEvent.theStripData), thePixelData(trackerEvent.thePixelData), thePhase2OTData(nullptr), + thePhase2OTVectorHits(nullptr), + thePhase2OTVectorHitsRej(nullptr), theOwner(false) { //std::cout << "Creatign non-owned MT @ " << this << " from @ " << & trackerEvent << " (strip data @ " << trackerEvent.theStripData << ")" << std::endl; if (stripClustersToSkip.refProd().id() != theStripData->handle().id()) { @@ -77,6 +83,8 @@ MeasurementTrackerEvent::MeasurementTrackerEvent( theStripData(nullptr), thePixelData(trackerEvent.thePixelData), thePhase2OTData(trackerEvent.thePhase2OTData), + thePhase2OTVectorHits(trackerEvent.thePhase2OTVectorHits), + thePhase2OTVectorHitsRej(trackerEvent.thePhase2OTVectorHitsRej), theOwner(false) { if (pixelClustersToSkip.refProd().id() != thePixelData->handle().id()) { edm::LogError("ProductIdMismatch") << "The pixel masking does not point to the proper collection of clusters: " diff --git a/RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h b/RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h index abe2be323bde5..893f7fdf1b2b9 100644 --- a/RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h +++ b/RecoTracker/TkDetLayers/interface/GeometricSearchTrackerBuilder.h @@ -18,7 +18,8 @@ class GeometricSearchTrackerBuilder { GeometricSearchTracker* build(const GeometricDet* theGeometricTracker, const TrackerGeometry* theGeomDetGeometry, - const TrackerTopology* tTopo) __attribute__((cold)); + const TrackerTopology* tTopo, + const bool usePhase2Stacks = false) __attribute__((cold)); }; #endif diff --git a/RecoTracker/TkDetLayers/src/GeometricSearchTrackerBuilder.cc b/RecoTracker/TkDetLayers/src/GeometricSearchTrackerBuilder.cc index 9be742143e173..51cd455a23717 100644 --- a/RecoTracker/TkDetLayers/src/GeometricSearchTrackerBuilder.cc +++ b/RecoTracker/TkDetLayers/src/GeometricSearchTrackerBuilder.cc @@ -21,7 +21,8 @@ using namespace std; GeometricSearchTracker* GeometricSearchTrackerBuilder::build(const GeometricDet* theGeometricTracker, const TrackerGeometry* theGeomDetGeometry, - const TrackerTopology* tTopo) { + const TrackerTopology* tTopo, + const bool usePhase2Stacks) { PixelBarrelLayerBuilder aPixelBarrelLayerBuilder; Phase2OTBarrelLayerBuilder aPhase2OTBarrelLayerBuilder; PixelForwardLayerBuilder aPixelForwardLayerBuilder; @@ -41,6 +42,7 @@ GeometricSearchTracker* GeometricSearchTrackerBuilder::build(const GeometricDet* vector thePosTIDLayers; vector theNegTECLayers; vector thePosTECLayers; + bool useBrothers = !usePhase2Stacks; vector theGeometricDetLayers = theGeometricTracker->components(); for (vector::const_iterator it = theGeometricDetLayers.begin(); @@ -93,10 +95,11 @@ GeometricSearchTracker* GeometricSearchTrackerBuilder::build(const GeometricDet* if ((*it)->type() == GeometricDet::OTPhase2Barrel) { vector theTOBGeometricDetLayers = (*it)->components(); + for (vector::const_iterator it2 = theTOBGeometricDetLayers.begin(); it2 != theTOBGeometricDetLayers.end(); it2++) { - theTOBLayers.push_back(aPhase2OTBarrelLayerBuilder.build(*it2, theGeomDetGeometry)); + theTOBLayers.push_back(aPhase2OTBarrelLayerBuilder.build(*it2, theGeomDetGeometry, useBrothers)); } } @@ -164,13 +167,15 @@ GeometricSearchTracker* GeometricSearchTrackerBuilder::build(const GeometricDet* if ((*it)->type() == GeometricDet::OTPhase2EndCap) { vector theTIDGeometricDetLayers = (*it)->components(); + + bool useBrothers = !usePhase2Stacks; for (vector::const_iterator it2 = theTIDGeometricDetLayers.begin(); it2 != theTIDGeometricDetLayers.end(); it2++) { if ((*it2)->positionBounds().z() < 0) - theNegTIDLayers.push_back(aPhase2EndcapLayerBuilder.build(*it2, theGeomDetGeometry, true)); + theNegTIDLayers.push_back(aPhase2EndcapLayerBuilder.build(*it2, theGeomDetGeometry, useBrothers)); if ((*it2)->positionBounds().z() > 0) - thePosTIDLayers.push_back(aPhase2EndcapLayerBuilder.build(*it2, theGeomDetGeometry, true)); + thePosTIDLayers.push_back(aPhase2EndcapLayerBuilder.build(*it2, theGeomDetGeometry, useBrothers)); } } diff --git a/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.cc b/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.cc index dd794dcab9b87..bd8beefb265dd 100644 --- a/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.cc +++ b/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.cc @@ -6,7 +6,7 @@ using namespace std; Phase2EndcapLayer* Phase2EndcapLayerBuilder::build(const GeometricDet* aPhase2EndcapLayer, const TrackerGeometry* theGeomDetGeometry, - const bool isOuterTracker) { + const bool useBrothers) { LogTrace("TkDetLayers") << "Phase2EndcapLayerBuilder::build"; vector theGeometricRings = aPhase2EndcapLayer->components(); LogTrace("TkDetLayers") << "theGeometricRings.size(): " << theGeometricRings.size(); @@ -18,8 +18,8 @@ Phase2EndcapLayer* Phase2EndcapLayerBuilder::build(const GeometricDet* aPhase2En it++) { // if we are in the phaseII OT, it will use the brothers to build pt modules // if we are in the phaseII pixel detector, it will not - thePhase2EndcapRings.push_back(myBuilder.build(*it, theGeomDetGeometry, isOuterTracker)); + thePhase2EndcapRings.push_back(myBuilder.build(*it, theGeomDetGeometry, useBrothers)); } - return new Phase2EndcapLayer(thePhase2EndcapRings, isOuterTracker); + return new Phase2EndcapLayer(thePhase2EndcapRings, useBrothers); } diff --git a/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.h b/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.h index 5a61c12b4dabe..b577b5b8c6380 100644 --- a/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.h +++ b/RecoTracker/TkDetLayers/src/Phase2EndcapLayerBuilder.h @@ -16,7 +16,7 @@ class Phase2EndcapLayerBuilder { Phase2EndcapLayerBuilder(){}; Phase2EndcapLayer* build(const GeometricDet* aPhase2EndcapLayer, const TrackerGeometry* theGeomDetGeometry, - const bool isOuterTracker) __attribute__((cold)); + const bool useBrothers) __attribute__((cold)); }; #pragma GCC visibility pop diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.cc b/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.cc index b950df541c2b9..81c43d72e3158 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.cc +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.cc @@ -7,7 +7,8 @@ using namespace std; using namespace edm; Phase2OTBarrelLayer* Phase2OTBarrelLayerBuilder::build(const GeometricDet* aPhase2OTBarrelLayer, - const TrackerGeometry* theGeomDetGeometry) { + const TrackerGeometry* theGeomDetGeometry, + const bool useBrothers) { // This builder is very similar to TOBLayer one. Most of the code should be put in a // common place. @@ -40,14 +41,17 @@ Phase2OTBarrelLayer* Phase2OTBarrelLayerBuilder::build(const GeometricDet* aPhas for (unsigned int index = 0; index != theGeometricDetRods.size(); index++) meanR += theGeometricDetRods[index]->positionBounds().perp(); if (!theGeometricDetRods.empty()) + meanR /= (double)theGeometricDetRods.size(); for (unsigned int index = 0; index != theGeometricDetRods.size(); index++) { if (theGeometricDetRods[index]->positionBounds().perp() < meanR) - theInnerRods.push_back(myPhase2OTBarrelRodBuilder.build(theGeometricDetRods[index], theGeomDetGeometry)); + theInnerRods.push_back( + myPhase2OTBarrelRodBuilder.build(theGeometricDetRods[index], theGeomDetGeometry, useBrothers)); if (theGeometricDetRods[index]->positionBounds().perp() > meanR) - theOuterRods.push_back(myPhase2OTBarrelRodBuilder.build(theGeometricDetRods[index], theGeomDetGeometry)); + theOuterRods.push_back( + myPhase2OTBarrelRodBuilder.build(theGeometricDetRods[index], theGeomDetGeometry, useBrothers)); } if (theGeometricDetRings.empty()) @@ -66,9 +70,9 @@ Phase2OTBarrelLayer* Phase2OTBarrelLayerBuilder::build(const GeometricDet* aPhas for (vector::const_iterator it = theGeometricDetRings.begin(); it != theGeometricDetRings.end(); it++) { if ((*it)->positionBounds().z() < centralZ) - theNegativeRings.push_back(myPhase2EndcapRingBuilder.build(*it, theGeomDetGeometry, true)); + theNegativeRings.push_back(myPhase2EndcapRingBuilder.build(*it, theGeomDetGeometry, useBrothers)); if ((*it)->positionBounds().z() > centralZ) - thePositiveRings.push_back(myPhase2EndcapRingBuilder.build(*it, theGeomDetGeometry, true)); + thePositiveRings.push_back(myPhase2EndcapRingBuilder.build(*it, theGeomDetGeometry, useBrothers)); } return new Phase2OTtiltedBarrelLayer(theInnerRods, theOuterRods, theNegativeRings, thePositiveRings); diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.h b/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.h index aa1ebc03af197..340193c50cc71 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.h +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelLayerBuilder.h @@ -14,8 +14,9 @@ class Phase2OTBarrelLayerBuilder { public: Phase2OTBarrelLayerBuilder(){}; - Phase2OTBarrelLayer* build(const GeometricDet* aPhase2OTBarrelLayer, const TrackerGeometry* theGeomDetGeometry) - __attribute__((cold)); + Phase2OTBarrelLayer* build(const GeometricDet* aPhase2OTBarrelLayer, + const TrackerGeometry* theGeomDetGeometry, + const bool useBrothers = true) __attribute__((cold)); }; #pragma GCC visibility pop diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.cc b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.cc index ba4a1bf7f5de2..b9cb3e4a22185 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.cc +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.cc @@ -26,8 +26,8 @@ namespace { Phase2OTBarrelRod::Phase2OTBarrelRod(vector& innerDets, vector& outerDets, - vector& innerDetBrothers, - vector& outerDetBrothers) + const vector& innerDetBrothers, + const vector& outerDetBrothers) : DetRod(true), theInnerDets(innerDets), theOuterDets(outerDets), @@ -61,6 +61,10 @@ Phase2OTBarrelRod::Phase2OTBarrelRod(vector& innerDets, << " , " << (**i).position().perp() << " , " << (**i).position().eta() << " , " << (**i).position().phi(); } + if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty()) + LogDebug("TkDetLayers") << "==== with stacks ====="; + if (!theInnerDetBrothers.empty() && !theOuterDetBrothers.empty()) + LogDebug("TkDetLayers") << "==== without stacks ====="; for (vector::const_iterator i = theOuterDets.begin(); i != theOuterDets.end(); i++) { LogDebug("TkDetLayers") << "outer Phase2OTBarrelRod's Det pos z,perp,eta,phi: " << (**i).position().z() << " , " @@ -198,6 +202,9 @@ bool Phase2OTBarrelRod::addClosest(const TrajectoryStateOnSurface& tsos, vector& brotherresult) const { const vector& sRod(subRod(crossing.subLayerIndex())); bool firstgroup = CompatibleDetToGroupAdder::add(*sRod[crossing.closestDetIndex()], tsos, prop, est, result); + if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty()) + return firstgroup; + // it assumes that the closestDetIndex is ok also for the brother detectors: the crossing is NOT recomputed const vector& sRodBrothers(subRodBrothers(crossing.subLayerIndex())); bool brothergroup = @@ -268,6 +275,8 @@ void Phase2OTBarrelRod::searchNeighbors(const TrajectoryStateOnSurface& tsos, break; if (!Adder::add(*sRod[idet], tsos, prop, est, result)) break; + if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty()) + break; // If the two above checks are passed also the brother module will be added with no further checks Adder::add(*sBrotherRod[idet], tsos, prop, est, brotherresult); } @@ -276,6 +285,8 @@ void Phase2OTBarrelRod::searchNeighbors(const TrajectoryStateOnSurface& tsos, break; if (!Adder::add(*sRod[idet], tsos, prop, est, result)) break; + if (theInnerDetBrothers.empty() && theOuterDetBrothers.empty()) + break; // If the two above checks are passed also the brother module will be added with no further checks Adder::add(*sBrotherRod[idet], tsos, prop, est, brotherresult); } diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.h b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.h index 3e0b77fcd3ea4..cf14dec6ac06d 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.h +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRod.h @@ -17,8 +17,9 @@ class Phase2OTBarrelRod final : public DetRod { Phase2OTBarrelRod(std::vector& innerDets, std::vector& outerDets, - std::vector& innerDetBrothers, - std::vector& outerDetBrothers) __attribute__((cold)); + const std::vector& innerDetBrothers = std::vector(), + const std::vector& outerDetBrothers = std::vector()) + __attribute__((cold)); ~Phase2OTBarrelRod() override __attribute__((cold)); // GeometricSearchDet interface diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.cc b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.cc index 786ffc5cd67a7..71d3286bdb129 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.cc +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.cc @@ -5,10 +5,11 @@ using namespace edm; using namespace std; Phase2OTBarrelRod* Phase2OTBarrelRodBuilder::build(const GeometricDet* thePhase2OTBarrelRod, - const TrackerGeometry* theGeomDetGeometry) { + const TrackerGeometry* theGeomDetGeometry, + const bool useBrothers) { vector allGeometricDets = thePhase2OTBarrelRod->components(); - vector compGeometricDets; LogDebug("TkDetLayers") << "Phase2OTBarrelRodBuilder with #Modules: " << allGeometricDets.size() << std::endl; + LogDebug("TkDetLayers") << " useBrothers: " << useBrothers << std::endl; vector innerGeomDets; vector outerGeomDets; @@ -16,47 +17,69 @@ Phase2OTBarrelRod* Phase2OTBarrelRodBuilder::build(const GeometricDet* thePhase2 vector outerGeomDetBrothers; double meanR = 0; - double meanRBrothers = 0; - for (vector::const_iterator it = allGeometricDets.begin(); it != allGeometricDets.end(); it++) { - compGeometricDets = (*it)->components(); - if (compGeometricDets.size() != 2) { - LogDebug("TkDetLayers") << " Stack not with two components but with " << compGeometricDets.size() << std::endl; - } else { - //LogTrace("TkDetLayers") << " compGeometricDets[0]->positionBounds().perp() " << compGeometricDets[0]->positionBounds().perp() << std::endl; - //LogTrace("TkDetLayers") << " compGeometricDets[1]->positionBounds().perp() " << compGeometricDets[1]->positionBounds().perp() << std::endl; - meanR = meanR + compGeometricDets[0]->positionBounds().perp(); - meanRBrothers = meanRBrothers + compGeometricDets[1]->positionBounds().perp(); + + if (!useBrothers) { + for (auto const& compGeometricDets : allGeometricDets) { + meanR = meanR + compGeometricDets->positionBounds().perp(); } - } - meanR = meanR / allGeometricDets.size(); - meanRBrothers = meanRBrothers / allGeometricDets.size(); - LogDebug("TkDetLayers") << " meanR Lower " << meanR << std::endl; - LogDebug("TkDetLayers") << " meanR Upper " << meanRBrothers << std::endl; + meanR = meanR / allGeometricDets.size(); + LogDebug("TkDetLayers") << " meanR Lower " << meanR << std::endl; + for (auto const& compGeometricDets : allGeometricDets) { + const GeomDet* theGeomDet = theGeomDetGeometry->idToDet(compGeometricDets->geographicalId()); + + if (compGeometricDets->positionBounds().perp() < meanR) + innerGeomDets.push_back(theGeomDet); - for (vector::iterator it = allGeometricDets.begin(); it != allGeometricDets.end(); it++) { - compGeometricDets = (*it)->components(); - const GeomDet* theGeomDet = theGeomDetGeometry->idToDet(compGeometricDets[0]->geographicalId()); - LogTrace("TkDetLayers") << " inserisco " << compGeometricDets[0]->geographicalId().rawId() << std::endl; + if (compGeometricDets->positionBounds().perp() > meanR) + outerGeomDets.push_back(theGeomDet); + } - if (compGeometricDets[0]->positionBounds().perp() < meanR) - innerGeomDets.push_back(theGeomDet); + LogDebug("TkDetLayers") << "innerGeomDets.size(): " << innerGeomDets.size(); + LogDebug("TkDetLayers") << "outerGeomDets.size(): " << outerGeomDets.size(); - if (compGeometricDets[0]->positionBounds().perp() > meanR) - outerGeomDets.push_back(theGeomDet); + } else { + vector compGeometricDets; - const GeomDet* theGeomDetBrother = theGeomDetGeometry->idToDet(compGeometricDets[1]->geographicalId()); - LogTrace("TkDetLayers") << " inserisco " << compGeometricDets[1]->geographicalId().rawId() << std::endl; - if (compGeometricDets[1]->positionBounds().perp() < meanRBrothers) - innerGeomDetBrothers.push_back(theGeomDetBrother); + double meanRBrothers = 0; + for (auto& it : allGeometricDets) { + compGeometricDets = it->components(); + if (compGeometricDets.size() != 2) { + LogDebug("TkDetLayers") << " Stack not with two components but with " << compGeometricDets.size() << std::endl; + } else { + meanR = meanR + compGeometricDets[0]->positionBounds().perp(); + meanRBrothers = meanRBrothers + compGeometricDets[1]->positionBounds().perp(); + } + } + meanR = meanR / allGeometricDets.size(); + meanRBrothers = meanRBrothers / allGeometricDets.size(); + LogDebug("TkDetLayers") << " meanR Lower " << meanR << std::endl; + LogDebug("TkDetLayers") << " meanR Upper " << meanRBrothers << std::endl; - if (compGeometricDets[1]->positionBounds().perp() > meanRBrothers) - outerGeomDetBrothers.push_back(theGeomDetBrother); - } + for (auto& it : allGeometricDets) { + compGeometricDets = it->components(); + const GeomDet* theGeomDet = theGeomDetGeometry->idToDet(compGeometricDets[0]->geographicalId()); + LogTrace("TkDetLayers") << " inserting " << compGeometricDets[0]->geographicalId().rawId() << std::endl; - LogDebug("TkDetLayers") << "innerGeomDets.size(): " << innerGeomDets.size(); - LogDebug("TkDetLayers") << "outerGeomDets.size(): " << outerGeomDets.size(); - LogDebug("TkDetLayers") << "innerGeomDetsBro.size(): " << innerGeomDetBrothers.size(); - LogDebug("TkDetLayers") << "outerGeomDetsBro.size(): " << outerGeomDetBrothers.size(); + if (compGeometricDets[0]->positionBounds().perp() < meanR) + innerGeomDets.push_back(theGeomDet); + + else + outerGeomDets.push_back(theGeomDet); + + const GeomDet* theGeomDetBrother = theGeomDetGeometry->idToDet(compGeometricDets[1]->geographicalId()); + LogTrace("TkDetLayers") << " inserting " << compGeometricDets[1]->geographicalId().rawId() << std::endl; + if (compGeometricDets[1]->positionBounds().perp() < meanRBrothers) + innerGeomDetBrothers.push_back(theGeomDetBrother); + + else + outerGeomDetBrothers.push_back(theGeomDetBrother); + } + + LogDebug("TkDetLayers") << "innerGeomDets.size(): " << innerGeomDets.size(); + LogDebug("TkDetLayers") << "outerGeomDets.size(): " << outerGeomDets.size(); + LogDebug("TkDetLayers") << "innerGeomDetsBro.size(): " << innerGeomDetBrothers.size(); + LogDebug("TkDetLayers") << "outerGeomDetsBro.size(): " << outerGeomDetBrothers.size(); + } return new Phase2OTBarrelRod(innerGeomDets, outerGeomDets, innerGeomDetBrothers, outerGeomDetBrothers); } diff --git a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.h b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.h index 641d04d4e7368..f51e0852451bb 100644 --- a/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.h +++ b/RecoTracker/TkDetLayers/src/Phase2OTBarrelRodBuilder.h @@ -14,8 +14,9 @@ class Phase2OTBarrelRodBuilder { public: Phase2OTBarrelRodBuilder(){}; - Phase2OTBarrelRod* build(const GeometricDet* thePhase2OTBarrelRod, const TrackerGeometry* theGeomDetGeometry) - __attribute__((cold)); + Phase2OTBarrelRod* build(const GeometricDet* thePhase2OTBarrelRod, + const TrackerGeometry* theGeomDetGeometry, + const bool useBrothers = true) __attribute__((cold)); }; #pragma GCC visibility pop diff --git a/RecoTracker/TkSeedGenerator/plugins/SeedingOTEDProducer.cc b/RecoTracker/TkSeedGenerator/plugins/SeedingOTEDProducer.cc new file mode 100644 index 0000000000000..ba28198988b29 --- /dev/null +++ b/RecoTracker/TkSeedGenerator/plugins/SeedingOTEDProducer.cc @@ -0,0 +1,372 @@ + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/Event.h" + +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" +#include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h" +#include "MagneticField/Engine/interface/MagneticField.h" + +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" +#include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h" + +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h" + +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +class TrajectoryStateUpdator; + +class SeedingOTEDProducer : public edm::stream::EDProducer<> { +public: + explicit SeedingOTEDProducer(const edm::ParameterSet&); + ~SeedingOTEDProducer() override; + void produce(edm::Event&, const edm::EventSetup&) override; + + static void fillDescriptions(edm::ConfigurationDescriptions&); + + TrajectorySeedCollection run(edm::Handle); + unsigned int checkLayer(unsigned int iidd); + std::vector collectVHsOnLayer(const edmNew::DetSetVector&, unsigned int); + void printVHsOnLayer(const edmNew::DetSetVector&, unsigned int); + const TrajectoryStateOnSurface buildInitialTSOS(const VectorHit*) const; + AlgebraicSymMatrix55 assign44To55(const AlgebraicSymMatrix44&) const; + std::pair propagateAndUpdate(const TrajectoryStateOnSurface initialTSOS, + const Propagator&, + const TrackingRecHit& hit) const; + float computeGlobalThetaError(const VectorHit* vh, const double sigmaZ_beamSpot) const; + float computeInverseMomentumError(const VectorHit* vh, + const float globalTheta, + const double sigmaZ_beamSpot, + const double transverseMomentum) const; + + TrajectorySeed createSeed(const TrajectoryStateOnSurface& tsos, + const edm::OwnVector& container, + const DetId& id, + const Propagator& prop) const; + + struct isInvalid { + bool operator()(const TrajectoryMeasurement& measurement) { + return ((measurement.recHit() == nullptr) || !(measurement.recHit()->isValid()) || + !((measurement).updatedState().isValid())); + } + }; + +private: + edm::EDGetTokenT vhProducerToken_; + const TrackerTopology* tkTopo_; + const MeasurementTracker* measurementTracker_; + std::unique_ptr layerMeasurements_; + const MeasurementEstimator* estimator_; + const Propagator* propagator_; + const MagneticField* magField_; + const TrajectoryStateUpdator* updator_; + const edm::EDGetTokenT tkMeasEventToken_; + edm::EDGetTokenT beamSpotToken_; + const reco::BeamSpot* beamSpot_; + std::string updatorName_; + + edm::ESGetToken topoToken_; + edm::ESGetToken propagatorToken_; + edm::ESGetToken magFieldToken_; + edm::ESGetToken updatorToken_; + edm::ESGetToken measurementTrackerToken_; + edm::ESGetToken estToken_; + edm::EDPutTokenT putToken_; +}; + +SeedingOTEDProducer::SeedingOTEDProducer(edm::ParameterSet const& conf) + : tkMeasEventToken_(consumes(conf.getParameter("trackerEvent"))), + topoToken_(esConsumes()), + propagatorToken_(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))), + magFieldToken_(esConsumes()), + updatorToken_(esConsumes(edm::ESInputTag("", "KFUpdator"))), + measurementTrackerToken_(esConsumes()), + estToken_(esConsumes(edm::ESInputTag("", "Chi2"))) { + vhProducerToken_ = consumes(edm::InputTag(conf.getParameter("src"))); + beamSpotToken_ = consumes(conf.getParameter("beamSpotLabel")); + updatorName_ = conf.getParameter("updator"); + putToken_ = produces(); +} + +SeedingOTEDProducer::~SeedingOTEDProducer() {} + +void SeedingOTEDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("src", edm::InputTag("siPhase2VectorHits", "accepted")); + desc.add("trackerEvent", edm::InputTag("MeasurementTrackerEvent")); + desc.add("beamSpotLabel", edm::InputTag("offlineBeamSpot")); + desc.add("updator", std::string("KFUpdator")); + descriptions.add("SeedingOTEDProducer", desc); +} + +void SeedingOTEDProducer::produce(edm::Event& event, const edm::EventSetup& es) { + tkTopo_ = &es.getData(topoToken_); + + edm::ESHandle measurementTrackerHandle = es.getHandle(measurementTrackerToken_); + measurementTracker_ = measurementTrackerHandle.product(); + + auto const& measurementTrackerEvent = event.get(tkMeasEventToken_); + + layerMeasurements_ = std::make_unique(*measurementTrackerHandle, measurementTrackerEvent); + + estimator_ = &es.getData(estToken_); + + propagator_ = &es.getData(propagatorToken_); + + magField_ = &es.getData(magFieldToken_); + + updator_ = &es.getData(updatorToken_); + + beamSpot_ = &event.get(beamSpotToken_); + + // Get the vector hits + auto vhs = event.getHandle(vhProducerToken_); + + event.emplace(putToken_, run(vhs)); +} + +TrajectorySeedCollection SeedingOTEDProducer::run(edm::Handle vHs) { + TrajectorySeedCollection result; + //check if all the first three layers have VHs + std::vector vhSeedsL1 = collectVHsOnLayer(*vHs.product(), 1); + std::vector vhSeedsL2 = collectVHsOnLayer(*vHs.product(), 2); + std::vector vhSeedsL3 = collectVHsOnLayer(*vHs.product(), 3); + if (vhSeedsL1.empty() || vhSeedsL2.empty() || vhSeedsL3.empty()) { + return result; + } + + //seeds are built in the L3 of the OT + const BarrelDetLayer* barrelOTLayer2 = measurementTracker_->geometricSearchTracker()->tobLayers().at(1); + + //the search propag directiondepend on the sign of signZ*signPz, while the building is always the contrary + Propagator* searchingPropagator = &*propagator_->clone(); + Propagator* buildingPropagator = &*propagator_->clone(); + buildingPropagator->setPropagationDirection(alongMomentum); + + for (const auto& hitL3 : vhSeedsL3) { + //building a tsos out of a VectorHit + const TrajectoryStateOnSurface initialTSOS = buildInitialTSOS(hitL3); + float signZ = copysign(1.0, initialTSOS.globalPosition().z()); + float signPz = copysign(1.0, initialTSOS.globalMomentum().z()); + + //set the direction of the propagator + if (signZ * signPz > 0.0) + searchingPropagator->setPropagationDirection(oppositeToMomentum); + else + searchingPropagator->setPropagationDirection(alongMomentum); + + //find vHits in layer 2 + std::vector measurementsL2 = + layerMeasurements_->measurements(*barrelOTLayer2, initialTSOS, *searchingPropagator, *estimator_); + + std::vector::iterator measurementsL2end = + std::remove_if(measurementsL2.begin(), measurementsL2.end(), isInvalid()); + measurementsL2.erase(measurementsL2end, measurementsL2.end()); + + if (!measurementsL2.empty()) { + //not sure if building it everytime takes time/memory + const DetLayer* barrelOTLayer1 = measurementTracker_->geometricSearchTracker()->tobLayers().at(0); + + for (const auto& mL2 : measurementsL2) { + const TrackingRecHit* hitL2 = mL2.recHit().get(); + + //propagate to the L2 and update the TSOS + std::pair updatedTSOS = + propagateAndUpdate(initialTSOS, *searchingPropagator, *hitL2); + if (!updatedTSOS.first) + continue; + + //searching possible VHs in L1 + std::vector measurementsL1 = + layerMeasurements_->measurements(*barrelOTLayer1, updatedTSOS.second, *searchingPropagator, *estimator_); + std::vector::iterator measurementsL1end = + std::remove_if(measurementsL1.begin(), measurementsL1.end(), isInvalid()); + measurementsL1.erase(measurementsL1end, measurementsL1.end()); + + if (!measurementsL1.empty()) { + for (const auto& mL1 : measurementsL1) { + const TrackingRecHit* hitL1 = mL1.recHit().get(); + + //propagate to the L1 and update the TSOS + std::pair updatedTSOSL1 = + propagateAndUpdate(updatedTSOS.second, *searchingPropagator, *hitL1); + if (!updatedTSOSL1.first) + continue; + + //building trajectory inside-out + if (searchingPropagator->propagationDirection() == alongMomentum) { + buildingPropagator->setPropagationDirection(oppositeToMomentum); + } else { + buildingPropagator->setPropagationDirection(alongMomentum); + } + + updatedTSOSL1.second.rescaleError(100); + + TrajectoryStateOnSurface updatedTSOSL1_final = updator_->update(updatedTSOSL1.second, *hitL1); + if UNLIKELY (!updatedTSOSL1_final.isValid()) + continue; + std::pair updatedTSOSL2_final = + propagateAndUpdate(updatedTSOSL1_final, *buildingPropagator, *hitL2); + if (!updatedTSOSL2_final.first) + continue; + std::pair updatedTSOSL3_final = + propagateAndUpdate(updatedTSOSL2_final.second, *buildingPropagator, *hitL3); + if (!updatedTSOSL3_final.first) + continue; + + edm::OwnVector container; + container.push_back(hitL1->clone()); + container.push_back(hitL2->clone()); + container.push_back(hitL3->clone()); + + TrajectorySeed ts = + createSeed(updatedTSOSL3_final.second, container, hitL3->geographicalId(), *buildingPropagator); + result.push_back(ts); + } + } + } + } + } + result.shrink_to_fit(); + return result; +} + +unsigned int SeedingOTEDProducer::checkLayer(unsigned int iidd) { + StripSubdetector strip = StripSubdetector(iidd); + unsigned int subid = strip.subdetId(); + if (subid == StripSubdetector::TIB || subid == StripSubdetector::TOB) { + return tkTopo_->layer(iidd); + } + return 0; +} + +std::vector SeedingOTEDProducer::collectVHsOnLayer(const edmNew::DetSetVector& input, + unsigned int layerNumber) { + std::vector vHsOnLayer; + if (!input.empty()) { + for (const auto& DSViter : input) { + if (checkLayer(DSViter.id()) == layerNumber) { + for (const auto& vh : DSViter) { + vHsOnLayer.push_back(&vh); + } + } + } + } + + return vHsOnLayer; +} + +void SeedingOTEDProducer::printVHsOnLayer(const edmNew::DetSetVector& input, unsigned int layerNumber) { + if (!input.empty()) { + for (const auto& DSViter : input) { + for (const auto& vh : DSViter) { + if (checkLayer(DSViter.id()) == layerNumber) + LogTrace("SeedingOTEDProducer") << " VH in layer " << layerNumber << " >> " << vh; + } + } + } else { + LogTrace("SeedingOTEDProducer") << " No VHs in layer " << layerNumber << "."; + } +} + +const TrajectoryStateOnSurface SeedingOTEDProducer::buildInitialTSOS(const VectorHit* vHit) const { + // having fun with theta + Global3DVector gv(vHit->globalPosition().x(), vHit->globalPosition().y(), vHit->globalPosition().z()); + float theta = gv.theta(); + // gv transform to local (lv) + const Local3DVector lv(vHit->det()->surface().toLocal(gv)); + + //FIXME::charge is fine 1 every two times!! + GlobalPoint center(0.0, 0.0, 0.0); + int charge = 1; + // momentum is a signed quantity in this case + float mom = vHit->momentum(magField_->inTesla(center).z()); + float xPos = vHit->localPosition().x(); + float yPos = vHit->localPosition().y(); + float dx = vHit->localDirection().x(); + // for dy use second component of the lv renormalized to the z component + float dy = lv.y() / lv.z(); + + // Pz and Dz should have the same sign + float signPz = copysign(1.0, vHit->globalPosition().z()); + + LocalTrajectoryParameters ltpar2(charge / mom, dx, dy, xPos, yPos, signPz); + AlgebraicSymMatrix55 mat = assign44To55(vHit->covMatrix()); + // set the error on 1/p + mat[0][0] = pow(computeInverseMomentumError( + vHit, theta, beamSpot_->sigmaZ(), vHit->transverseMomentum(magField_->inTesla(center).z())), + 2); + + //building tsos + LocalTrajectoryError lterr(mat); + const TrajectoryStateOnSurface tsos(ltpar2, lterr, vHit->det()->surface(), magField_); + + return tsos; +} + +AlgebraicSymMatrix55 SeedingOTEDProducer::assign44To55(const AlgebraicSymMatrix44& mat44) const { + AlgebraicSymMatrix55 result; + for (int i = 1; i < 5; i++) { + for (int j = 1; j < 5; j++) { + result[i][j] = mat44[i - 1][j - 1]; + } + } + return result; +} + +std::pair SeedingOTEDProducer::propagateAndUpdate( + const TrajectoryStateOnSurface initialTSOS, const Propagator& prop, const TrackingRecHit& hit) const { + TrajectoryStateOnSurface propTSOS = prop.propagate(initialTSOS, hit.det()->surface()); + if UNLIKELY (!propTSOS.isValid()) + return std::make_pair(false, propTSOS); + TrajectoryStateOnSurface updatedTSOS = updator_->update(propTSOS, hit); + if UNLIKELY (!updatedTSOS.isValid()) + return std::make_pair(false, updatedTSOS); + return std::make_pair(true, updatedTSOS); +} + +float SeedingOTEDProducer::computeGlobalThetaError(const VectorHit* vh, const double sigmaZ_beamSpot) const { + double derivative = vh->globalPosition().perp() / vh->globalPosition().mag2(); + double derivative2 = pow(derivative, 2); + return pow(derivative2 * vh->lowerGlobalPosErr().czz() + derivative2 * pow(sigmaZ_beamSpot, 2), 0.5); +} + +float SeedingOTEDProducer::computeInverseMomentumError(const VectorHit* vh, + const float globalTheta, + const double sigmaZ_beamSpot, + const double transverseMomentum) const { + //for pT > 2GeV, 1/pT has sigma = 1/sqrt(12) + float varianceInverseTransvMomentum = 1. / 12; + float derivativeTheta2 = pow(cos(globalTheta) / transverseMomentum, 2); + float derivativeInverseTransvMomentum2 = pow(sin(globalTheta), 2); + float thetaError = computeGlobalThetaError(vh, sigmaZ_beamSpot); + return pow(derivativeTheta2 * pow(thetaError, 2) + derivativeInverseTransvMomentum2 * varianceInverseTransvMomentum, + 0.5); +} + +TrajectorySeed SeedingOTEDProducer::createSeed(const TrajectoryStateOnSurface& tsos, + const edm::OwnVector& container, + const DetId& id, + const Propagator& prop) const { + PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId()); + return TrajectorySeed(seedTSOS, container, prop.propagationDirection()); +} +DEFINE_FWK_MODULE(SeedingOTEDProducer); diff --git a/RecoTracker/TkSeedingLayers/src/HitExtractorSTRP.cc b/RecoTracker/TkSeedingLayers/src/HitExtractorSTRP.cc index 44aa1e9d10cbd..471c131a5f4c0 100644 --- a/RecoTracker/TkSeedingLayers/src/HitExtractorSTRP.cc +++ b/RecoTracker/TkSeedingLayers/src/HitExtractorSTRP.cc @@ -39,6 +39,7 @@ HitExtractorSTRP::HitExtractorSTRP(GeomDetEnumerators::SubDetector subdet, hasMatchedHits(false), hasRPhiHits(false), hasStereoHits(false), + hasVectorHits(false), hasRingSelector(false), hasSimpleRphiHitsCleaner(true) { minGoodCharge = iminGoodCharge; @@ -48,6 +49,7 @@ HitExtractorSTRP::HitExtractorSTRP(GeomDetEnumerators::SubDetector subdet, void HitExtractorSTRP::useSkipClusters_(const edm::InputTag& m, edm::ConsumesCollector& iC) { theSkipClusters = iC.consumes(m); + theSkipPhase2Clusters = iC.consumes(m); } void HitExtractorSTRP::useRingSelector(int minRing, int maxRing) { @@ -69,8 +71,8 @@ bool HitExtractorSTRP::skipThis( if (maskCluster && (stripClusterMask->mask(clus.key()))) return true; - if - UNLIKELY(minGoodCharge <= 0) return false; + if UNLIKELY (minGoodCharge <= 0) + return false; return siStripClusterTools::chargePerCM(id, *clus.cluster_strip()) <= minGoodCharge; } @@ -116,49 +118,74 @@ void HitExtractorSTRP::cleanedOfClusters(const TkTransientTrackingRecHitBuilder& HitExtractor::Hits& hits, bool matched, unsigned int cleanFrom) const { - LogDebug("HitExtractorPIX") << "getting: " << hits.size() << " in input."; - edm::Handle stripClusterMask; - if (maskCluster) - ev.getByToken(theSkipClusters, stripClusterMask); unsigned int skipped = 0; unsigned int projected = 0; - for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) { - assert(hits[iH]->isValid()); - auto id = hits[iH]->geographicalId(); - if (matched) { - bool replace; - ProjectedSiStripRecHit2D* replaceMe; - std::tie(replace, replaceMe) = skipThis(ttrhBuilder, *hits[iH], stripClusterMask); - if (replace) { - if (!replaceMe) { - LogDebug("HitExtractorSTRP") << "skipping a matched hit on :" << hits[iH]->geographicalId().rawId(); - skipped++; - } else - projected++; - hits[iH].reset(replaceMe); - if (replaceMe == nullptr) - assert(hits[iH].empty()); - else - assert(hits[iH].isOwn()); + if (hasMatchedHits || hasRPhiHits || hasStereoHits) { + LogTrace("HitExtractorSTRP") << "getting " << hits.size() << " strip hit in input."; + edm::Handle stripClusterMask; + if (maskCluster) + ev.getByToken(theSkipClusters, stripClusterMask); + for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) { + assert(hits[iH]->isValid()); + auto id = hits[iH]->geographicalId(); + if (matched) { + auto [replace, replaceMe] = skipThis(ttrhBuilder, *hits[iH], stripClusterMask); + if (replace) { + if (!replaceMe) { + LogTrace("HitExtractorSTRP") << "skipping a matched hit on :" << hits[iH]->geographicalId().rawId(); + skipped++; + } else { + projected++; + } + hits[iH].reset(replaceMe); + if (replaceMe == nullptr) + assert(hits[iH].empty()); + else + assert(hits[iH].isOwn()); + } + } else if (skipThis(id, hits[iH]->firstClusterRef(), stripClusterMask)) { + LogTrace("HitExtractorSTRP") << "skipping a hit on :" << hits[iH]->geographicalId().rawId() << " key: "; + skipped++; + hits[iH].reset(); + } + } + } + if (hasVectorHits) { + LogTrace("HitExtractorSTRP") << "getting " << hits.size() << " vector hit in input."; + edm::Handle ph2ClusterMask; + if (maskCluster) + ev.getByToken(theSkipPhase2Clusters, ph2ClusterMask); + for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) { + LogTrace("HitExtractorSTRP") << "analizing hit on :" << hits[iH]->geographicalId().rawId(); + assert(hits[iH]->isValid()); + const VectorHit& vhit = dynamic_cast(*hits[iH]); + LogTrace("HitExtractorSTRP") << " key lower: " << vhit.lowerClusterRef().key() + << " and key upper: " << vhit.upperClusterRef().key(); + LogTrace("HitExtractorSTRP") << " key lower: " << hits[iH]->firstClusterRef().key(); + + //FIXME:: introduce a "projected" version later? + if (maskCluster && + (ph2ClusterMask->mask(vhit.lowerClusterRef().key()) || ph2ClusterMask->mask(vhit.upperClusterRef().key()))) { + LogTrace("HitExtractorSTRP") << "skipping a vector hit on :" << hits[iH]->geographicalId().rawId() + << " key lower: " << vhit.lowerClusterRef().key() + << " and key upper: " << vhit.upperClusterRef().key(); + skipped++; + hits[iH].reset(); } - } else if (skipThis(id, hits[iH]->firstClusterRef(), stripClusterMask)) { - LogDebug("HitExtractorSTRP") << "skipping a hit on :" << hits[iH]->geographicalId().rawId() << " key: "; - skipped++; - hits[iH].reset(); } } + // remove empty elements... auto last = std::remove_if(hits.begin() + cleanFrom, hits.end(), [](HitPointer const& p) { return p.empty(); }); hits.resize(last - hits.begin()); - // std::cout << "HitExtractorSTRP " <<"skipped :"< matchedHits; ev.getByToken(theMatchedHits, matchedHits); @@ -208,12 +236,24 @@ HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder if (skipClusters) cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); } + if (hasVectorHits) { + LogError("HitExtractorSTRP") << "TIB is not supposed to be in Phase2 TRK detector configuration. What follows " + "have never been checked before! "; + auto const& vectorHits = ev.get(theVectorHits); + if (skipClusters) + cleanFrom = result.size(); + range2SeedingHits(vectorHits, result, tTopo->tibDetIdLayerComparator(theIdLayer)); + if (skipClusters) + cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); + } + } // // TID // else if (theLayerSubDet == GeomDetEnumerators::TID) { + LogTrace("HitExtractorSTRP") << "Getting hits into the TID"; if (hasMatchedHits) { edm::Handle matchedHits; ev.getByToken(theMatchedHits, matchedHits); @@ -271,11 +311,32 @@ HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder if (skipClusters) cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); } + if (hasVectorHits) { + LogTrace("HitExtractorSTRP") << "Getting vector hits for IdLayer " << theIdLayer; + auto const& vectorHits = ev.get(theVectorHits); + //FIXME: check the skipClusters with VHits + if (skipClusters) + cleanFrom = result.size(); + auto getter = tTopo->tidDetIdWheelComparator(static_cast(theSide), theIdLayer); + VectorHitCollection::Range range = vectorHits.equal_range(getter.first, getter.second); + for (VectorHitCollection::const_iterator it = range.first; it != range.second; ++it) { + int ring = tTopo->tidRing(it->detId()); + if (!ringRange(ring)) + continue; + for (VectorHitCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) { + result.emplace_back(*hit); + } + } + LogTrace("HitExtractorSTRP") << "result size value:" << result.size(); + if (skipClusters) + cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); + } } // // TOB // else if (theLayerSubDet == GeomDetEnumerators::TOB) { + LogTrace("HitExtractorSTRP") << "Getting hits into the TOB"; if (hasMatchedHits) { edm::Handle matchedHits; ev.getByToken(theMatchedHits, matchedHits); @@ -325,12 +386,25 @@ HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder if (skipClusters) cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); } + if (hasVectorHits) { + LogTrace("HitExtractorSTRP") << "Getting vector hits for IdLayer " << theIdLayer; + edm::Handle vectorHits; + ev.getByToken(theVectorHits, vectorHits); + //FIXME: check the skipClusters with VHits + if (skipClusters) + cleanFrom = result.size(); + range2SeedingHits(*vectorHits, result, tTopo->tobDetIdLayerComparator(theIdLayer)); + if (skipClusters) + cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); + } + } // // TEC // else if (theLayerSubDet == GeomDetEnumerators::TEC) { + LogTrace("HitExtractorSTRP") << "Getting hits into the TEC"; if (hasMatchedHits) { edm::Handle matchedHits; ev.getByToken(theMatchedHits, matchedHits); @@ -388,9 +462,28 @@ HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder if (skipClusters) cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); } + if (hasVectorHits) { + LogError("HitExtractorSTRP") << "TEC is not supposed to be in Phase2 TRK detector configuration. What follows " + "have never been checked before! "; + edm::Handle vectorHits; + ev.getByToken(theVectorHits, vectorHits); + if (skipClusters) + cleanFrom = result.size(); + auto getter = tTopo->tidDetIdWheelComparator(static_cast(theSide), theIdLayer); + VectorHitCollection::Range range = vectorHits->equal_range(getter.first, getter.second); + for (VectorHitCollection::const_iterator it = range.first; it != range.second; ++it) { + int ring = tTopo->tidRing(it->detId()); + if (!ringRange(ring)) + continue; + for (VectorHitCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) { + result.emplace_back(*hit); + } + } + if (skipClusters) + cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom); + } } - LogDebug("HitExtractorSTRP") << " giving: " << result.size() << " out"; - // std::cout << "HitExtractorSTRP "<<" giving: "< #include @@ -44,6 +45,11 @@ namespace ctfseeding { hasStereoHits = true; theStereoHits = iC.consumes(m); } + + void useVectorHits(const edm::InputTag& m, edm::ConsumesCollector& iC) { + hasVectorHits = true; + theVectorHits = iC.consumes(m); + } void useRingSelector(int minRing, int maxRing); void useSimpleRphiHitsCleaner(bool use) { hasSimpleRphiHitsCleaner = use; } @@ -72,6 +78,7 @@ namespace ctfseeding { bool ringRange(int ring) const; typedef edm::ContainerMask > SkipClustersCollection; + typedef edm::ContainerMask SkipPhase2ClustersCollection; void useSkipClusters_(const edm::InputTag& m, edm::ConsumesCollector& iC) override; private: @@ -81,12 +88,15 @@ namespace ctfseeding { double minAbsZ; int theMinRing, theMaxRing; edm::EDGetTokenT theSkipClusters; + edm::EDGetTokenT theSkipPhase2Clusters; edm::EDGetTokenT theMatchedHits; edm::EDGetTokenT theRPhiHits; edm::EDGetTokenT theStereoHits; + edm::EDGetTokenT theVectorHits; bool hasMatchedHits; bool hasRPhiHits; bool hasStereoHits; + bool hasVectorHits; bool hasRingSelector; bool hasSimpleRphiHitsCleaner; bool failProjection; diff --git a/RecoTracker/TkSeedingLayers/src/SeedingLayerSetsBuilder.cc b/RecoTracker/TkSeedingLayers/src/SeedingLayerSetsBuilder.cc index 96b4ed9599a93..4bd23772b1761 100644 --- a/RecoTracker/TkSeedingLayers/src/SeedingLayerSetsBuilder.cc +++ b/RecoTracker/TkSeedingLayers/src/SeedingLayerSetsBuilder.cc @@ -133,6 +133,9 @@ SeedingLayerSetsBuilder::LayerSpec::LayerSpec(unsigned short index, if (cfgLayer.exists("stereoRecHits")) { extr->useStereoHits(cfgLayer.getParameter("stereoRecHits"), iC); } + if (cfgLayer.exists("vectorRecHits")) { + extr->useVectorHits(cfgLayer.getParameter("vectorRecHits"), iC); + } if (cfgLayer.exists("useRingSlector") && cfgLayer.getParameter("useRingSlector")) { extr->useRingSelector(cfgLayer.getParameter("minRing"), cfgLayer.getParameter("maxRing")); } diff --git a/RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h b/RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h index d12a0bf069028..9fdecece36a57 100644 --- a/RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h +++ b/RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h @@ -33,6 +33,7 @@ class TkClonerImpl final : public TkCloner { TrajectoryStateOnSurface const& tsos) const override; std::unique_ptr operator()(Phase2TrackerRecHit1D const& hit, TrajectoryStateOnSurface const& tsos) const override; + std::unique_ptr operator()(VectorHit const& hit, TrajectoryStateOnSurface const& tsos) const override; using TkCloner::makeShared; TrackingRecHit::ConstRecHitPointer makeShared(SiPixelRecHit const& hit, @@ -47,7 +48,8 @@ class TkClonerImpl final : public TkCloner { TrajectoryStateOnSurface const& tsos) const override; TrackingRecHit::ConstRecHitPointer makeShared(Phase2TrackerRecHit1D const& hit, TrajectoryStateOnSurface const& tsos) const override; - + TrackingRecHit::ConstRecHitPointer makeShared(VectorHit const& hit, + TrajectoryStateOnSurface const& tsos) const override; // project either mono or stero hit... std::unique_ptr project(SiStripMatchedRecHit2D const& hit, bool mono, diff --git a/RecoTracker/TransientTrackingRecHit/src/TkClonerImpl.cc b/RecoTracker/TransientTrackingRecHit/src/TkClonerImpl.cc index b8c2b24f19888..8c396591f7a96 100644 --- a/RecoTracker/TransientTrackingRecHit/src/TkClonerImpl.cc +++ b/RecoTracker/TransientTrackingRecHit/src/TkClonerImpl.cc @@ -6,6 +6,7 @@ #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" @@ -56,6 +57,10 @@ std::unique_ptr TkClonerImpl::operator()(Phase2TrackerRec return std::make_unique(params.first, params.second, *hit.det(), hit.cluster()); } +std::unique_ptr TkClonerImpl::operator()(VectorHit const& hit, TrajectoryStateOnSurface const& tsos) const { + return std::make_unique(hit); +} + TrackingRecHit::ConstRecHitPointer TkClonerImpl::makeShared(SiPixelRecHit const& hit, TrajectoryStateOnSurface const& tsos) const { // std::cout << "cloning " << typeid(hit).name() << std::endl; @@ -92,6 +97,11 @@ TrackingRecHit::ConstRecHitPointer TkClonerImpl::makeShared(Phase2TrackerRecHit1 return std::make_unique(params.first, params.second, *hit.det(), hit.cluster()); } +TrackingRecHit::ConstRecHitPointer TkClonerImpl::makeShared(VectorHit const& hit, + TrajectoryStateOnSurface const& tsos) const { + return std::make_shared(hit); +} + namespace { #undef RecoTracker_TransientTrackingRecHit_TSiStripMatchedRecHit_RefitProj #undef RecoTracker_TransientTrackingRecHit_TSiStripMatchedRecHit_RefitLGL diff --git a/SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h b/SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h index 3a17d4ab3e874..76a56e9f0ebd4 100644 --- a/SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h +++ b/SimTracker/TrackAssociation/interface/trackHitsToClusterRefs.h @@ -8,6 +8,7 @@ #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/VectorHit.h" #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h" namespace track_associator { @@ -56,6 +57,13 @@ namespace track_associator { edm::LogError("TrackAssociator") << ">>> RecHit does not have an associated cluster!" << " file: " << __FILE__ << " line: " << __LINE__; returnValue.push_back(ph2Hit->omniClusterRef()); + } else if (tid == typeid(VectorHit)) { + const VectorHit *vectorHit = dynamic_cast(rhit); + if (!vectorHit->cluster().isNonnull()) + edm::LogError("TrackAssociator") << ">>> RecHit does not have an associated cluster!" + << " file: " << __FILE__ << " line: " << __LINE__; + returnValue.push_back(vectorHit->firstClusterRef()); + } else { auto const &thit = static_cast(*rhit); if (thit.isProjected()) { diff --git a/SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.cc b/SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.cc index f22a5c69e4352..2b1f01b89c6f4 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.cc +++ b/SimTracker/TrackAssociatorProducers/plugins/QuickTrackAssociatorByHitsImpl.cc @@ -801,6 +801,5 @@ double QuickTrackAssociatorByHitsImpl::weightedNumberOfTrackClusters(iter begin, } } LogTrace("QuickTrackAssociatorByHitsImpl") << " total weighted clusters: " << weightedClusters; - return weightedClusters; }