diff --git a/Configuration/EventContent/python/EventContent_cff.py b/Configuration/EventContent/python/EventContent_cff.py index 7f5f99daf1c64..67912ad685215 100644 --- a/Configuration/EventContent/python/EventContent_cff.py +++ b/Configuration/EventContent/python/EventContent_cff.py @@ -900,3 +900,10 @@ def _addOutputCommands(mod, newCommands): _addOutputCommands(FEVTEventContent,RecoLocalFastTimeFEVT) _addOutputCommands(RECOSIMEventContent,RecoLocalFastTimeRECO) _addOutputCommands(AODSIMEventContent,RecoLocalFastTimeAOD) + +from RecoMTD.Configuration.RecoMTD_EventContent_cff import RecoMTDFEVT, RecoMTDRECO, RecoMTDAOD +_addOutputCommands(FEVTDEBUGEventContent,RecoMTDFEVT) +_addOutputCommands(FEVTDEBUGHLTEventContent,RecoMTDFEVT) +_addOutputCommands(FEVTEventContent,RecoMTDFEVT) +_addOutputCommands(RECOSIMEventContent,RecoMTDRECO) +_addOutputCommands(AODSIMEventContent,RecoMTDAOD) diff --git a/Configuration/StandardSequences/python/Reconstruction_cff.py b/Configuration/StandardSequences/python/Reconstruction_cff.py index 822da0e74a98c..bc63e0f0240cf 100644 --- a/Configuration/StandardSequences/python/Reconstruction_cff.py +++ b/Configuration/StandardSequences/python/Reconstruction_cff.py @@ -5,6 +5,7 @@ from RecoLocalMuon.Configuration.RecoLocalMuon_cff import * from RecoLocalCalo.Configuration.RecoLocalCalo_cff import * from RecoLocalFastTime.Configuration.RecoLocalFastTime_cff import * +from RecoMTD.Configuration.RecoMTD_cff import * from RecoTracker.Configuration.RecoTracker_cff import * from RecoParticleFlow.PFClusterProducer.particleFlowCluster_cff import * from TrackingTools.Configuration.TrackingTools_cff import * @@ -124,6 +125,10 @@ _run3_globalreco = globalreco.copyAndExclude([CastorFullReco]) run3_common.toReplaceWith(globalreco, _run3_globalreco) +_phase2_timing_layer_globalreco = _run3_globalreco.copy() +_phase2_timing_layer_globalreco += fastTimingGlobalReco +phase2_timing_layer.toReplaceWith(globalreco,_phase2_timing_layer_globalreco) + _fastSim_globalreco = globalreco.copyAndExclude([CastorFullReco,muoncosmicreco]) # insert the few tracking modules to be run after mixing back in the globalreco sequence _fastSim_globalreco.insert(0,newCombinedSeeds+trackExtrapolator+caloTowerForTrk+firstStepPrimaryVerticesUnsorted+ak4CaloJetsForTrk+initialStepTrackRefsForJets+firstStepPrimaryVertices) diff --git a/DQM/TrackingMonitor/python/trackingRecoMaterialAnalyzer_cfi.py b/DQM/TrackingMonitor/python/trackingRecoMaterialAnalyzer_cfi.py index 3a1c329068264..31831d8618352 100644 --- a/DQM/TrackingMonitor/python/trackingRecoMaterialAnalyzer_cfi.py +++ b/DQM/TrackingMonitor/python/trackingRecoMaterialAnalyzer_cfi.py @@ -11,6 +11,7 @@ TrackerRecHitBuilder = cms.string('WithAngleAndTemplate'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite'), diff --git a/DataFormats/FTLRecHit/interface/FTLCluster.h b/DataFormats/FTLRecHit/interface/FTLCluster.h new file mode 100644 index 0000000000000..eabf978f8f46c --- /dev/null +++ b/DataFormats/FTLRecHit/interface/FTLCluster.h @@ -0,0 +1,272 @@ +#ifndef DataFormats_FTLRecHit_FTLCluster_h +#define DataFormats_FTLRecHit_FTLCluster_h + +/** \class FTLCluster + * + * based on SiPixelCluster + * + * \author Paolo Meridiani + */ + +#include +#include +#include +#include +#include +#include + +#include "DataFormats/DetId/interface/DetId.h" + +class FTLCluster { +public: + + typedef DetId key_type; + + class FTLHit { + public: + constexpr FTLHit() : x_(0), y_(0), energy_(0), time_(0), time_error_(0) {} + constexpr FTLHit(uint16_t hit_x, uint16_t hit_y, float hit_energy, float hit_time, float hit_time_error) : + x_(hit_x), y_(hit_y), energy_(hit_energy), time_(hit_time), time_error_(hit_time_error) {} + constexpr uint16_t x() { return x_; } + constexpr uint16_t y() { return y_; } + constexpr uint16_t energy() { return energy_; } + constexpr uint16_t time() { return time_; } + constexpr uint16_t time_error() { return time_error_; } + private: + uint16_t x_; //row + uint16_t y_; //col + float energy_; + float time_; + float time_error_; + }; + + //--- Integer shift in x and y directions. + class Shift { + public: + constexpr Shift( int dx, int dy) : dx_(dx), dy_(dy) {} + constexpr Shift() : dx_(0), dy_(0) {} + constexpr int dx() const { return dx_;} + constexpr int dy() const { return dy_;} + private: + int dx_; + int dy_; + }; + + //--- Position of a FTL Hit + class FTLHitPos { + public: + constexpr FTLHitPos() : row_(0), col_(0) {} + constexpr FTLHitPos(int row, int col) : row_(row) , col_(col) {} + constexpr int row() const { return row_;} + constexpr int col() const { return col_;} + constexpr FTLHitPos operator+( const Shift& shift) const { + return FTLHitPos( row() + shift.dx(), col() + shift.dy()); + } + private: + int row_; + int col_; + }; + + static constexpr unsigned int MAXSPAN=255; + static constexpr unsigned int MAXPOS=2047; + + /** Construct from a range of digis that form a cluster and from + * a DetID. The range is assumed to be non-empty. + */ + FTLCluster() {} + + FTLCluster(DetId id, unsigned int isize, float const * energys, float const* times, float const* time_errors, + uint16_t const * xpos, uint16_t const * ypos, + uint16_t const xmin, uint16_t const ymin) : + theid(id), theHitOffset(2*isize), theHitENERGY(energys,energys+isize), theHitTIME(times,times+isize), theHitTIME_ERROR(time_errors,time_errors+isize) { + uint16_t maxCol = 0; + uint16_t maxRow = 0; + int maxHit=-1; + float maxEnergy=-99999; + for (unsigned int i=0; i!=isize; ++i) { + uint16_t xoffset = xpos[i]-xmin; + uint16_t yoffset = ypos[i]-ymin; + theHitOffset[i*2] = std::min(uint16_t(MAXSPAN),xoffset); + theHitOffset[i*2+1] = std::min(uint16_t(MAXSPAN),yoffset); + if (xoffset > maxRow) maxRow = xoffset; + if (yoffset > maxCol) maxCol = yoffset; + if (theHitENERGY[i]>maxEnergy) + { + maxHit=i; + maxEnergy=theHitENERGY[i]; + } + } + packRow(xmin,maxRow); + packCol(ymin,maxCol); + + if (maxHit>=0) + seed_=std::min(uint8_t(MAXSPAN),uint8_t(maxHit)); + } + + // linear average position (barycenter) + inline float x() const { + auto x_pos=[this](unsigned int i) { return this->theHitOffset[i*2] + minHitRow() + 0.5f; }; + return weighted_mean(this->theHitENERGY,x_pos); + } + + inline float y() const { + auto y_pos=[this](unsigned int i) { return this->theHitOffset[i*2+1] + minHitCol() + 0.5f; }; + return weighted_mean(this->theHitENERGY,y_pos); + } + + inline float time() const { + auto t=[this](unsigned int i) { return this->theHitTIME[i]; }; + return weighted_mean(this->theHitENERGY,t); + } + + inline float timeError() const { + auto t_err=[this](unsigned int i) { return this->theHitTIME_ERROR[i]; }; + return weighted_mean_error(this->theHitENERGY,t_err); + } + + // Return number of hits. + inline int size() const { return theHitENERGY.size();} + + // Return cluster dimension in the x direction. + inline int sizeX() const { return rowSpan() +1;} + + // Return cluster dimension in the y direction. + inline int sizeY() const { return colSpan() +1;} + + + inline float energy() const { + float qm = 0; + int isize = theHitENERGY.size(); + for (int i=0; i & hitOffset() const { return theHitOffset;} + const std::vector & hitENERGY() const { return theHitENERGY;} + const std::vector & hitTIME() const { return theHitTIME;} + const std::vector & hitTIME_ERROR() const { return theHitTIME_ERROR;} + + // infinite faster than above... + FTLHit hit(int i) const { + return FTLHit(minHitRow() + theHitOffset[i*2], + minHitCol() + theHitOffset[i*2+1], + theHitENERGY[i], + theHitTIME[i], + theHitTIME_ERROR[i] + ); + + } + + FTLHit seed() const { + return hit(seed_); + } + + int colSpan() const {return theHitColSpan; } + + int rowSpan() const { return theHitRowSpan; } + + const DetId& id() const { return theid; } + const DetId& detid() const { return id(); } + + bool overflowCol() const { return overflow_(theHitColSpan); } + + bool overflowRow() const { return overflow_(theHitRowSpan); } + + bool overflow() const { return overflowCol() || overflowRow(); } + + void packCol(uint16_t ymin, uint16_t yspan) { + theMinHitCol = ymin; + theHitColSpan = std::min(yspan, uint16_t(MAXSPAN)); + } + void packRow(uint16_t xmin, uint16_t xspan) { + theMinHitRow = xmin; + theHitRowSpan = std::min(xspan, uint16_t(MAXSPAN)); + } + + void setClusterErrorX( float errx ) { err_x = errx; } + void setClusterErrorY( float erry ) { err_y = erry; } + void setClusterErrorTime( float errtime ) { err_time = errtime; } + float getClusterErrorX() const { return err_x; } + float getClusterErrorY() const { return err_y; } + float getClusterErrorTime() const { return err_time; } + +private: + + DetId theid; + + std::vector theHitOffset; + std::vector theHitENERGY; + std::vector theHitTIME; + std::vector theHitTIME_ERROR; + + uint16_t theMinHitRow=MAXPOS; // Minimum hit index in the x direction (low edge). + uint16_t theMinHitCol=MAXPOS; // Minimum hit index in the y direction (left edge). + uint8_t theHitRowSpan=0; // Span hit index in the x direction (low edge). + uint8_t theHitColSpan=0; // Span hit index in the y direction (left edge). + + float err_x=-99999.9f; + float err_y=-99999.9f; + float err_time=-99999.9f; + + uint8_t seed_; + + float weighted_sum(const std::vector& weights, const std::function& sumFunc, const std::function& outFunc) const + { + float tot=0; + float sumW=0; + for (unsigned int i=0; i& weights, const std::function& value) const + { + auto sumFunc=[weights,value](unsigned int i) { return weights[i]*value(i); } ; + auto outFunc=[](float x,float y) { if (y>0) return (float)x/y; else return -999.f; }; + return weighted_sum(weights,sumFunc,outFunc); + } + + float weighted_mean_error(const std::vector& weights, const std::function& err) const + { + auto sumFunc=[weights,err](unsigned int i) { return weights[i]*weights[i]*err(i)*err(i); } ; + auto outFunc=[](float x,float y) { if (y>0) return (float)sqrt(x)/y; else return -999.f; }; + return weighted_sum(weights,sumFunc,outFunc); + } + + static int overflow_(uint16_t span) { return span==uint16_t(MAXSPAN);} + +}; + + +// Comparison operators (needed by DetSetVector & SortedCollection ) +inline bool operator<( const FTLCluster& one, const FTLCluster& other) { + if(one.detid() == other.detid()) { + if ( one.minHitRow() < other.minHitRow() ) { + return true; + } else if ( one.minHitRow() > other.minHitRow() ) { + return false; + } else if ( one.minHitCol() < other.minHitCol() ) { + return true; + } else { + return false; + } + } + return one.detid() < other.detid(); +} + +inline bool operator<( const FTLCluster& one, const uint32_t& detid) { + return one.detid() < detid;} + +inline bool operator<( const uint32_t& detid, const FTLCluster& other) { + return detid < other.detid();} + +#endif diff --git a/DataFormats/FTLRecHit/interface/FTLClusterCollections.h b/DataFormats/FTLRecHit/interface/FTLClusterCollections.h new file mode 100644 index 0000000000000..1c277139d8fe4 --- /dev/null +++ b/DataFormats/FTLRecHit/interface/FTLClusterCollections.h @@ -0,0 +1,17 @@ +#ifndef DataFormats_FTLRecHit_FTLClusterCollections_H +#define DataFormats_FTLRecHit_FTLClusterCollections_H + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/Common/interface/DetSetRefVector.h" +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/Common/interface/RefVector.h" + +#include "DataFormats/FTLRecHit/interface/FTLCluster.h" + +typedef edmNew::DetSetVector FTLClusterCollection; +typedef edm::Ref FTLClusterRef; +typedef edm::DetSetRefVector FTLClusterRefs; +typedef edm::RefProd FTLClustersRef; + +#endif diff --git a/DataFormats/FTLRecHit/src/FTLCluster.cc b/DataFormats/FTLRecHit/src/FTLCluster.cc new file mode 100644 index 0000000000000..e03a27dfe3de0 --- /dev/null +++ b/DataFormats/FTLRecHit/src/FTLCluster.cc @@ -0,0 +1 @@ +#include "DataFormats/FTLRecHit/interface/FTLCluster.h" diff --git a/DataFormats/FTLRecHit/src/classes.h b/DataFormats/FTLRecHit/src/classes.h index fda7aecfe7591..086366d5be36b 100644 --- a/DataFormats/FTLRecHit/src/classes.h +++ b/DataFormats/FTLRecHit/src/classes.h @@ -1,8 +1,11 @@ #include "DataFormats/FTLRecHit/interface/FTLUncalibratedRecHit.h" #include "DataFormats/FTLRecHit/interface/FTLRecHit.h" #include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h" -#include "DataFormats/FTLRecHit/interface/FTLSeverityLevel.h" + +#include "DataFormats/FTLRecHit/interface/FTLCluster.h" +#include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h" #include "DataFormats/FTLRecHit/interface/FTLTrackingRecHit.h" +#include "DataFormats/FTLRecHit/interface/FTLSeverityLevel.h" #include "DataFormats/Common/interface/RefProd.h" #include "DataFormats/Common/interface/Wrapper.h" #include "DataFormats/Common/interface/RefToBase.h" @@ -30,12 +33,24 @@ namespace DataFormats_FTLRecHit { FTLRecHitRef _FTLRHitRef; FTLRecHitRefs _FTLRHitRefs; FTLRecHitsRef _FTLRHitsRef; - + + FTLCluster _aCluster; + std::vector _FTLClusterVect; + FTLClusterCollection _theFTLCsc; + + edm::Wrapper< FTLClusterCollection > _FTLClusterProd; + FTLClusterRef _FTLClusterRef; + FTLClusterRefs _FTLClusterRefs; + FTLClustersRef _FTLClustersRef; + + std::vector,FTLCluster,edmNew::DetSetVector::FindForDetSetVector> > dsvr_v; + edmNew::DetSetVector,FTLCluster,edmNew::DetSetVector::FindForDetSetVector> > dsvr; + edm::Wrapper,FTLCluster,edmNew::DetSetVector::FindForDetSetVector> > > dsvr_w; + FTLTrackingRecHitFromHit _FTLTTrackingRecHitFromRef; FTLTrackingRecHitCollection _FTLTTrackingRecHitCollection; edm::Wrapper _FTLTTrackingRecHitProd; - - + }; } @@ -54,5 +69,13 @@ namespace DataFormats_FTLRecHit { edm::Wrapper< std::vector > > dummy31; edm::Wrapper< edm::DetSetVector > dummy41; edm::Wrapper< std::vector< std::vector < edm::DetSet > > > dummy51; + + edm::Wrapper< FTLCluster > dummy02; + edm::Wrapper< std::vector > dummy12; + edm::Wrapper< edm::DetSet > dummy22; + edm::Wrapper< std::vector > > dummy32; + edm::Wrapper< edm::DetSetVector > dummy42; + edm::Wrapper< std::vector< std::vector < edm::DetSet > > > dummy52; + }; } diff --git a/DataFormats/FTLRecHit/src/classes_def.xml b/DataFormats/FTLRecHit/src/classes_def.xml index 7308af710c1dc..cb10113fc9e3c 100644 --- a/DataFormats/FTLRecHit/src/classes_def.xml +++ b/DataFormats/FTLRecHit/src/classes_def.xml @@ -1,7 +1,7 @@ - + @@ -20,6 +20,16 @@ + + + + + + + + + + @@ -30,7 +40,14 @@ + + + + + + + diff --git a/DataFormats/TrackReco/interface/TrackBase.h b/DataFormats/TrackReco/interface/TrackBase.h index cecaae7859c26..77617eaad0737 100644 --- a/DataFormats/TrackReco/interface/TrackBase.h +++ b/DataFormats/TrackReco/interface/TrackBase.h @@ -171,7 +171,7 @@ class TrackBase const Vector &momentum, int charge, const CovarianceMatrix &cov, TrackAlgorithm = undefAlgorithm, TrackQuality quality = undefQuality, signed char nloops = 0, uint8_t stopReason = 0, - float t0 = 0.f, float beta = 0., + float t0 = 0.f, float beta = 0.f, float covt0t0 = -1.f, float covbetabeta = -1.f); /// virtual destructor diff --git a/DataFormats/TrackerRecHit2D/BuildFile.xml b/DataFormats/TrackerRecHit2D/BuildFile.xml index 0d1b1b6dff0d3..b669fdbfdcc49 100644 --- a/DataFormats/TrackerRecHit2D/BuildFile.xml +++ b/DataFormats/TrackerRecHit2D/BuildFile.xml @@ -1,5 +1,6 @@ + diff --git a/DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h b/DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h new file mode 100644 index 0000000000000..35edd27084331 --- /dev/null +++ b/DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h @@ -0,0 +1,41 @@ +#ifndef DataFormats_TrackerRecHit2D_MTDTrackingRecHit_h +#define DataFormats_TrackerRecHit2D_MTDTrackingRecHit_h + +/// A 2D TrackerRecHit with time and time error information + +#include +#include "DataFormats/TrackerRecHit2D/interface/TrackerSingleRecHit.h" +#include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h" + +class MTDTrackingRecHit : public TrackerSingleRecHit { + public: + + MTDTrackingRecHit() : TrackerSingleRecHit() {} + + MTDTrackingRecHit(const LocalPoint& p, const LocalError& e, + const GeomDet& idet, const FTLClusterRef& objref) : + TrackerSingleRecHit(p, e, idet, trackerHitRTTI::mipTiming, objref) + {} + + MTDTrackingRecHit* clone() const override { return new MTDTrackingRecHit(*this); } + + // things to specialize from BaseTrackerRecHit + bool isPhase2() const final { return true; } + void getKfComponents(KfComponentsHolder& holder) const final; + + int dimension() const final { return 2; } + + //specific timing stuff + float energy() const { return omniCluster().mtdCluster().energy(); } + float time() const { return omniCluster().mtdCluster().time(); } + float timeError() const { return omniCluster().mtdCluster().timeError(); } +}; + + +// Instantiations and specializations for FTLRecHitRef and reco::CaloClusterPtr +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/OwnVector.h" +typedef edmNew::DetSetVector MTDTrackingDetSetVector; +typedef edm::OwnVector MTDTrackingOwnVector; + +#endif diff --git a/DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h b/DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h index 644e1752815ca..fa477a9fd73b5 100644 --- a/DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h +++ b/DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h @@ -7,28 +7,32 @@ #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h" #include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h" class OmniClusterRef { - static const unsigned int kInvalid = 0x80000000; // bit 31 on - static const unsigned int kIsStrip = 0x20000000; // bit 29 on + static const unsigned int kInvalid = 0x80000000; // bit 31 on + static const unsigned int kIsStrip = 0x20000000; // bit 29 on // FIXME:: need to check when introducing phase2 pixel static const unsigned int kIsPhase2 = 0x40000000; // bit 30 on + static const unsigned int kIsTiming = 0x10000000; // bit 28 on static const unsigned int kIsRegional = 0x60000000; // bit 30 and 29 on (will become fastsim???) - - static const unsigned int indexMask = 0xFFFFFF; - static const unsigned int subClusMask = 0xF000000; + + static const unsigned int indexMask = 0xFFFFFF; + static const unsigned int subClusMask = 0xF000000; static const unsigned int subClusShift = 24; public: typedef edm::Ref,SiPixelCluster > ClusterPixelRef; typedef edm::Ref,SiStripCluster > ClusterStripRef; typedef edm::Ref, Phase2TrackerCluster1D> Phase2Cluster1DRef; + typedef edm::Ref ClusterMTDRef; OmniClusterRef() : me(edm::RefCore(),kInvalid) {} explicit OmniClusterRef(ClusterPixelRef const & ref, unsigned int subClus=0) : me(ref.refCore(), (ref.isNonnull() ? ref.key() | (subClus< TrackerSingleRecHit(const LocalPoint& p, const LocalError& e, GeomDet const & idet, trackerHitRTTI::RTTI rt, @@ -64,6 +66,10 @@ class TrackerSingleRecHit : public BaseTrackerRecHit { return cluster_.cluster_phase2OT(); } + ClusterMTDRef cluster_mtd() const { + return cluster_.cluster_mtd(); + } + SiStripCluster const & stripCluster() const { return cluster_.stripCluster(); } @@ -76,10 +82,15 @@ class TrackerSingleRecHit : public BaseTrackerRecHit { return cluster_.phase2OTCluster(); } + FTLCluster const & mtdCluster() const { + return cluster_.mtdCluster(); + } + // void setClusterRef(const & OmniClusterRef ref) { cluster_ =ref;} void setClusterPixelRef(ClusterPixelRef const & ref) { cluster_ = OmniClusterRef(ref); } void setClusterStripRef(ClusterStripRef const & ref) { cluster_ = OmniClusterRef(ref); } void setClusterPhase2Ref(ClusterPhase2Ref const & ref) { cluster_ = OmniClusterRef(ref); } + void setClusterMTDRef(ClusterMTDRef const & ref) { cluster_ = OmniClusterRef(ref); } bool sharesInput( const TrackingRecHit* other, SharedInputType what) const final; diff --git a/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h b/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h index 21fd379c68c0c..025235d5000aa 100644 --- a/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h +++ b/DataFormats/TrackerRecHit2D/interface/trackerHitRTTI.h @@ -7,7 +7,7 @@ namespace trackerHitRTTI { // tracking hit can be : single (si1D, si2D, pix), projected, matched or multi enum RTTI { undef=0, single=1, projStereo=2, projMono=3, match=4, multi=5, fastSingle=6, fastProjStereo=7,fastProjMono=8,fastMatch=9, - notFromCluster=10}; + notFromCluster=10, mipTiming=11}; inline RTTI rtti(TrackingRecHit const & hit) { return RTTI(hit.getRTTI());} inline bool isUndef(TrackingRecHit const & hit) { return rtti(hit)==undef;} inline bool isNotFromCluster(TrackingRecHit const & hit) { return rtti(hit)==notFromCluster;} @@ -21,6 +21,7 @@ namespace trackerHitRTTI { inline bool isFromDet(TrackingRecHit const & hit) { return (rtti(hit)>0) & (rtti(hit)<6) ;} 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 isTiming(TrackingRecHit const & hit) { return rtti(hit)==mipTiming; } inline unsigned int projId(TrackingRecHit const & hit) { return hit.rawId()+int(rtti(hit))-1;} } diff --git a/DataFormats/TrackerRecHit2D/src/MTDTrackingRecHit.cc b/DataFormats/TrackerRecHit2D/src/MTDTrackingRecHit.cc new file mode 100644 index 0000000000000..46814b6e32e6d --- /dev/null +++ b/DataFormats/TrackerRecHit2D/src/MTDTrackingRecHit.cc @@ -0,0 +1,5 @@ +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" + +void MTDTrackingRecHit::getKfComponents(KfComponentsHolder& holder) const { + getKfComponents2D(holder); +} diff --git a/DataFormats/TrackerRecHit2D/src/classes.h b/DataFormats/TrackerRecHit2D/src/classes.h index 5444a8657c5a9..3b4f47e1a27b1 100644 --- a/DataFormats/TrackerRecHit2D/src/classes.h +++ b/DataFormats/TrackerRecHit2D/src/classes.h @@ -27,6 +27,7 @@ #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h" #include "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" #include namespace DataFormats_TrackerRecHit2D { @@ -36,12 +37,14 @@ namespace DataFormats_TrackerRecHit2D { SiStripRecHit1D a11; SiStripMatchedRecHit2D a2; SiPixelRecHit b1; + MTDTrackingRecHit mtd1; edm::ClonePolicy a4; edm::ClonePolicy a44; edm::ClonePolicy a5; edm::ClonePolicy b2; edm::ClonePolicy e2; + edm::ClonePolicy e10; edm::OwnVector > a6; @@ -63,6 +66,7 @@ namespace DataFormats_TrackerRecHit2D { edm::ClonePolicy > e3; edm::OwnVector >::const_iterator it10; + MTDTrackingOwnVector::const_iterator it11; edm::OwnVector ovbtrh; edm::Wrapper> wovbtrh; @@ -104,6 +108,15 @@ namespace DataFormats_TrackerRecHit2D { edm::ClonePolicy >, edm::ClonePolicy >::id_iterator itpix; + edm::Wrapper< edm::RangeMap >, + edm::ClonePolicy > > mtdRecHitCollectionWrapper; + edm::RangeMap >, + edm::ClonePolicy >::id_iterator mtdpix; + edm::Ref >,edm::ClonePolicy >,SiStripRecHit2D,edm::refhelper::FindUsingAdvance >,edm::ClonePolicy >,SiStripRecHit2D> > refRangeMapDetIdOwnVectorSiStripRecHit2D; edm::RefVector >,edm::ClonePolicy >,SiStripRecHit2D,edm::refhelper::FindUsingAdvance >,edm::ClonePolicy >,SiStripRecHit2D> > refVectorRangeMapDetIdOwnVectorSiStripRecHit2D; @@ -115,6 +128,7 @@ namespace DataFormats_TrackerRecHit2D { edm::Wrapper > wdstvDummy11; edm::Wrapper > wdstvDummy2; edm::Wrapper > wdstvDummy3; + edm::Wrapper wdstvDummy4; edm::Wrapper clusterRemovalInfo; diff --git a/DataFormats/TrackerRecHit2D/src/classes_def.xml b/DataFormats/TrackerRecHit2D/src/classes_def.xml index c4daab655997e..291b20d729255 100644 --- a/DataFormats/TrackerRecHit2D/src/classes_def.xml +++ b/DataFormats/TrackerRecHit2D/src/classes_def.xml @@ -176,4 +176,14 @@ + + + + + + + + + + diff --git a/Geometry/GlobalTrackingGeometryBuilder/plugins/GlobalTrackingGeometryESProducer.cc b/Geometry/GlobalTrackingGeometryBuilder/plugins/GlobalTrackingGeometryESProducer.cc index 570c2d987e955..4529fcd058b2d 100644 --- a/Geometry/GlobalTrackingGeometryBuilder/plugins/GlobalTrackingGeometryESProducer.cc +++ b/Geometry/GlobalTrackingGeometryBuilder/plugins/GlobalTrackingGeometryESProducer.cc @@ -36,73 +36,72 @@ GlobalTrackingGeometryESProducer::produce(const GlobalTrackingGeometryRecord& re GEMGeometry const* gem = nullptr; ME0Geometry const* me0 = nullptr; - try { - edm::ESHandle tkH; - record.getRecord().get(tkH); + edm::ESHandle tkH; + if( auto tkRecord = record.tryToGetRecord() ) { + tkRecord->get(tkH); if(tkH.isValid()) { tk = tkH.product(); } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No Tracker geometry is available."; } - } catch (edm::eventsetup::NoRecordException& e){ + } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No TrackerDigiGeometryRecord is available."; } - try { - edm::ESHandle mtdH; - record.getRecord().get(mtdH); - if(mtdH.isValid()) { + edm::ESHandle mtdH; + if( auto mtdRecord = record.tryToGetRecord() ) { + mtdRecord->get(mtdH); + if( mtdH.isValid() ) { mtd = mtdH.product(); } else { - LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No MTD geometry is available."; + LogInfo("GeometryGlobalTrackingGeometryBuilder") << "No MTD geometry is available."; } - } catch (edm::eventsetup::NoRecordException& e){ - LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No MTDDigiGeometryRecord is available."; + } else { + LogInfo("GeometryGlobalTrackingGeometryBuilder") << "No MTDDigiGeometryRecord is available."; } - - - try { - edm::ESHandle dtH; - record.getRecord().get(dtH); + + edm::ESHandle dtH; + edm::ESHandle cscH; + edm::ESHandle rpcH; + edm::ESHandle gemH; + edm::ESHandle me0H; + if( auto muonRecord = record.tryToGetRecord() ) { + muonRecord->get(dtH); if(dtH.isValid()) { dt = dtH.product(); } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No DT geometry is available."; } - - edm::ESHandle cscH; - record.getRecord().get(cscH); + + muonRecord->get(cscH); if(cscH.isValid()) { csc = cscH.product(); } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No CSC geometry is available."; - } + } - edm::ESHandle rpcH; - record.getRecord().get(rpcH); + muonRecord->get(rpcH); if(rpcH.isValid()) { rpc = rpcH.product(); } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No RPC geometry is available."; } - - edm::ESHandle gemH; - record.getRecord().get(gemH); + + muonRecord->get(gemH); if(gemH.isValid()) { gem = gemH.product(); } else { LogInfo("GeometryGlobalTrackingGeometryBuilder") << "No GEM geometry is available."; } - - edm::ESHandle me0H; - record.getRecord().get(me0H); + + muonRecord->get(me0H); if(me0H.isValid()) { me0 = me0H.product(); } else { LogInfo("GeometryGlobalTrackingGeometryBuilder") << "No ME0 geometry is available."; } - } catch (edm::eventsetup::NoRecordException& e){ + } else { LogWarning("GeometryGlobalTrackingGeometryBuilder") << "No MuonGeometryRecord is available."; } diff --git a/Geometry/GlobalTrackingGeometryBuilder/test/BuildFile.xml b/Geometry/GlobalTrackingGeometryBuilder/test/BuildFile.xml index 254eb302d50b7..a2eb4e977a135 100644 --- a/Geometry/GlobalTrackingGeometryBuilder/test/BuildFile.xml +++ b/Geometry/GlobalTrackingGeometryBuilder/test/BuildFile.xml @@ -1,5 +1,6 @@ + diff --git a/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_EventContent_cff.py b/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_EventContent_cff.py index 79e261dea2db2..00b42da701603 100644 --- a/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_EventContent_cff.py +++ b/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_EventContent_cff.py @@ -7,16 +7,35 @@ 'keep *_ftlRecHits_*_*', 'keep *_mtdUncalibratedRecHits_*_*', 'keep *_mtdRecHits_*_*', + 'keep *_mtdClusters_*_*', + 'keep *_mtdTrackingRecHits_*_*' ) ) #RECO content RecoLocalFastTimeRECO = cms.PSet( outputCommands = cms.untracked.vstring( - 'keep *_ftlRecHits_*_*', + 'keep *_ftlRecHits_*_*', + 'keep *_mtdRecHits_*_*', + 'keep *_mtdClusters_*_*', ) ) #AOD content RecoLocalFastTimeAOD = cms.PSet( - outputCommands = cms.untracked.vstring( + outputCommands = cms.untracked.vstring( + 'keep *_mtdClusters_*_*', ) ) + +from Configuration.Eras.Modifier_phase2_timing_layer_tile_cff import phase2_timing_layer_tile +from Configuration.Eras.Modifier_phase2_timing_layer_bar_cff import phase2_timing_layer_bar + +(phase2_timing_layer_tile | phase2_timing_layer_bar).toModify( + RecoLocalFastTimeFEVT.outputCommands, + func=lambda outputCommands: outputCommands.extend(['drop *_ftlUncalibratedRecHits_*_*', + 'drop *_ftlRecHits_*_*']) +) + +(phase2_timing_layer_tile | phase2_timing_layer_bar).toModify( + RecoLocalFastTimeRECO.outputCommands, + func=lambda outputCommands: outputCommands.extend(['drop *_ftlRecHits_*_*']) +) diff --git a/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_cff.py b/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_cff.py index 5535ad828af3c..d24b7d4d50701 100644 --- a/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_cff.py +++ b/RecoLocalFastTime/Configuration/python/RecoLocalFastTime_cff.py @@ -8,8 +8,12 @@ from RecoLocalFastTime.FTLRecProducers.mtdUncalibratedRecHits_cfi import mtdUncalibratedRecHits from RecoLocalFastTime.FTLRecProducers.mtdRecHits_cfi import mtdRecHits +from RecoLocalFastTime.FTLRecProducers.mtdTrackingRecHits_cfi import mtdTrackingRecHits +from RecoLocalFastTime.FTLClusterizer.mtdClusters_cfi import mtdClusters -_phase2_timing_layer_fastTimingLocalReco = cms.Sequence(mtdUncalibratedRecHits*mtdRecHits) +from RecoLocalFastTime.FTLClusterizer.MTDCPEESProducers_cff import * + +_phase2_timing_layer_fastTimingLocalReco = cms.Sequence(mtdUncalibratedRecHits*mtdRecHits*mtdClusters*mtdTrackingRecHits) from Configuration.Eras.Modifier_phase2_timing_layer_tile_cff import phase2_timing_layer_tile from Configuration.Eras.Modifier_phase2_timing_layer_bar_cff import phase2_timing_layer_bar diff --git a/RecoLocalFastTime/FTLClusterizer/BuildFile.xml b/RecoLocalFastTime/FTLClusterizer/BuildFile.xml new file mode 100644 index 0000000000000..31d4276b46cf8 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/BuildFile.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h new file mode 100644 index 0000000000000..fefb11548e6ef --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h @@ -0,0 +1,143 @@ +#ifndef RecoLocalTracker_FTLClusterizer_MTDArrayBuffer_H +#define RecoLocalTracker_FTLClusterizer_MTDArrayBuffer_H + +//---------------------------------------------------------------------------- +//! \class MTDArrayBuffer +//! \brief Class to store ADC counts and times during clustering. +//! +//---------------------------------------------------------------------------- + +// We use FTLHitPos which is an inner class of FTLCluster: +#include "DataFormats/FTLRecHit/interface/FTLCluster.h" + +#include +#include + +class MTDArrayBuffer +{ + public: + typedef unsigned int uint; + + inline MTDArrayBuffer( uint rows, uint cols); + inline MTDArrayBuffer( ){} + + inline void setSize( uint rows, uint cols); + + inline float energy( uint row, uint col) const; + inline float energy( const FTLCluster::FTLHitPos&) const; + inline float time( uint row, uint col) const; + inline float time( const FTLCluster::FTLHitPos&) const; + inline float time_error( uint row, uint col) const; + inline float time_error( const FTLCluster::FTLHitPos&) const; + + inline uint rows() const { return nrows;} + inline uint columns() const { return ncols;} + + inline bool inside(uint row, uint col) const; + + inline void clear(uint row, uint col) + { + set_energy( row, col, 0.); + set_time( row, col, 0.); + set_time_error( row, col, 0.); + } + inline void clear(const FTLCluster::FTLHitPos& pos) + { + clear(pos.row(),pos.col()); + } + + inline void set( uint row, uint col, float energy, float time, float time_error); + inline void set( const FTLCluster::FTLHitPos&, float energy, float time, float time_error); + + inline void set_energy( uint row, uint col, float energy); + inline void set_energy( const FTLCluster::FTLHitPos&, float energy); + inline void add_energy( uint row, uint col, float energy); + + inline void set_time( uint row, uint col, float time); + inline void set_time( const FTLCluster::FTLHitPos&, float time); + + inline void set_time_error( uint row, uint col, float time_error); + inline void set_time_error( const FTLCluster::FTLHitPos&, float time_error); + + uint size() const { return hitEnergy_vec.size();} + + /// Definition of indexing within the buffer. + uint index( uint row, uint col) const {return col*nrows+row;} + uint index( const FTLCluster::FTLHitPos& pix) const { return index(pix.row(), pix.col()); } + + private: + std::vector hitEnergy_vec; + std::vector hitTime_vec; + std::vector hitTimeError_vec; + uint nrows; + uint ncols; +}; + +MTDArrayBuffer::MTDArrayBuffer( uint rows, uint cols) + : hitEnergy_vec(rows*cols,0), hitTime_vec(rows*cols,0), hitTimeError_vec(rows*cols,0), nrows(rows), ncols(cols) {} + +void MTDArrayBuffer::setSize( uint rows, uint cols) { + hitEnergy_vec.resize(rows*cols,0); + hitTime_vec.resize(rows*cols,0); + hitTimeError_vec.resize(rows*cols,0); + nrows = rows; + ncols = cols; +} + +bool MTDArrayBuffer::inside(uint row, uint col) const +{ + return ( row < nrows && col < ncols); +} + +float MTDArrayBuffer::energy(uint row, uint col) const { return hitEnergy_vec[index(row,col)];} +float MTDArrayBuffer::energy(const FTLCluster::FTLHitPos& pix) const {return hitEnergy_vec[index(pix)];} + +float MTDArrayBuffer::time(uint row, uint col) const { return hitTime_vec[index(row,col)];} +float MTDArrayBuffer::time(const FTLCluster::FTLHitPos& pix) const {return hitTime_vec[index(pix)];} + +float MTDArrayBuffer::time_error(uint row, uint col) const { return hitTimeError_vec[index(row,col)];} +float MTDArrayBuffer::time_error(const FTLCluster::FTLHitPos& pix) const {return hitTimeError_vec[index(pix)];} + +void MTDArrayBuffer::set( uint row, uint col, float energy, float time, float time_error) +{ + hitEnergy_vec[index(row,col)] = energy; + hitTime_vec[index(row,col)] = time; + hitTimeError_vec[index(row,col)] = time_error; +} +void MTDArrayBuffer::set( const FTLCluster::FTLHitPos& pix, float energy, float time, float time_error) +{ + set( pix.row(), pix.col(), energy, time, time_error); +} + +void MTDArrayBuffer::set_energy( uint row, uint col, float energy) +{ + hitEnergy_vec[index(row,col)] = energy; +} +void MTDArrayBuffer::set_energy( const FTLCluster::FTLHitPos& pix, float energy) +{ + hitEnergy_vec[index(pix)] = energy; +} +void MTDArrayBuffer::add_energy( uint row, uint col, float energy) +{ + hitEnergy_vec[index(row,col)] += energy; +} + +void MTDArrayBuffer::set_time( uint row, uint col, float time) +{ + hitTime_vec[index(row,col)] = time; +} +void MTDArrayBuffer::set_time( const FTLCluster::FTLHitPos& pix, float time) +{ + hitTime_vec[index(pix)] = time; +} + +void MTDArrayBuffer::set_time_error( uint row, uint col, float time_error) +{ + hitTimeError_vec[index(row,col)] = time_error; +} +void MTDArrayBuffer::set_time_error( const FTLCluster::FTLHitPos& pix, float time_error) +{ + hitTimeError_vec[index(pix)] = time_error; +} + +#endif diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h new file mode 100644 index 0000000000000..22f5abc5f4108 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h @@ -0,0 +1,127 @@ +#ifndef RecoLocalFastTime_FTLClusterizer_MTDCPEBase_H +#define RecoLocalFastTime_FTLClusterizer_MTDCPEBase_H 1 + +//----------------------------------------------------------------------------- +// \class MTDCPEBase +//----------------------------------------------------------------------------- + +#include +#include +#include +#include "TMath.h" + +#include "MTDClusterParameterEstimator.h" + +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetUnit.h" +#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" +#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h" +#include "DataFormats/GeometrySurface/interface/GloballyPositioned.h" + + +#include + +#include + + +class MTDCPEBase : public MTDClusterParameterEstimator +{ +public: + struct DetParam + { + DetParam() {} + const MTDGeomDetUnit * theDet; + const ProxyMTDTopology * theTopol; + const RectangularMTDTopology * theRecTopol; + + GeomDetType::SubDetector thePart; + Local3DPoint theOrigin; + float theThickness; + float thePitchX; + float thePitchY; + }; + + struct ClusterParam + { + ClusterParam(const FTLCluster & cl) : theCluster(&cl) {} + + virtual ~ClusterParam() = default; + + const FTLCluster * theCluster; + }; + +public: + MTDCPEBase(edm::ParameterSet const& conf, const MTDGeometry& geom); + + + inline ReturnType getParameters(const FTLCluster & cl, + const GeomDetUnit & det ) const override + { + + DetParam const & dp = detParam(det); + ClusterParam cp(cl); + setTheClu( dp, cp ); + auto tuple = std::make_tuple( + localPosition(dp, cp), + localError(dp, cp), + clusterTime(dp, cp), + clusterTimeError(dp, cp) + ); + return tuple; + } + + //-------------------------------------------------------------------------- + // In principle we could use the track too to get an angle if needed + //-------------------------------------------------------------------------- + inline ReturnType getParameters(const FTLCluster & cl, + const GeomDetUnit & det, + const LocalTrajectoryParameters & ltp ) const override + { + return getParameters(cl,det); + } + + + +private: + //-------------------------------------------------------------------------- + // This is where the action happens. + //-------------------------------------------------------------------------- + virtual LocalPoint localPosition(DetParam const & dp, ClusterParam & cp) const; + virtual LocalError localError (DetParam const & dp, ClusterParam & cp) const; + virtual TimeValue clusterTime(DetParam const & dp, ClusterParam & cp) const; + virtual TimeValueError clusterTimeError(DetParam const & dp, ClusterParam & cp) const; + + void fillDetParams(); + +protected: + //--------------------------------------------------------------------------- + // Data members + //--------------------------------------------------------------------------- + + //--- Global quantities + const MTDGeometry & geom_; // geometry + +protected: + + void setTheClu( DetParam const & dp, ClusterParam & cp ) const ; + + //--------------------------------------------------------------------------- + // Cluster-level services. + //--------------------------------------------------------------------------- + + DetParam const & detParam(const GeomDetUnit & det) const; + + using DetParams=std::vector; + DetParams m_DetParams; + +}; + +#endif + + diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h new file mode 100644 index 0000000000000..aec8395bd3af0 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h @@ -0,0 +1,63 @@ +#ifndef RecoLocalFastTime_MTDCluster_Parameter_Estimator_H +#define RecoLocalFastTime_MTDCluster_Parameter_Estimator_H + +#include "DataFormats/GeometrySurface/interface/LocalError.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" + +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" + +#include "DataFormats/FTLRecHit/interface/FTLCluster.h" + +#include + +class MTDClusterParameterEstimator +{ + public: + + virtual ~MTDClusterParameterEstimator(){} + + typedef std::pair LocalValues; + typedef std::vector VLocalValues; + + typedef float TimeValue; + typedef float TimeValueError; + + using ReturnType = std::tuple; + + // here just to implement it in the clients; + // to be properly implemented in the sub-classes in order to make them thread-safe + + virtual ReturnType getParameters(const FTLCluster & cl, + const GeomDetUnit & det) const =0; + + virtual ReturnType getParameters(const FTLCluster & cl, + const GeomDetUnit & det, + const LocalTrajectoryParameters & ltp ) const =0; + + virtual ReturnType getParameters(const FTLCluster & cl, + const GeomDetUnit & det, + const TrajectoryStateOnSurface& tsos ) const { + return getParameters(cl,det,tsos.localParameters()); + } + + virtual VLocalValues localParametersV(const FTLCluster& cluster, const GeomDetUnit& gd) const { + VLocalValues vlp; + ReturnType tuple = getParameters(cluster, gd); + vlp.emplace_back(std::get<0>(tuple), std::get<1>(tuple)); + return vlp; + } + virtual VLocalValues localParametersV(const FTLCluster& cluster, const GeomDetUnit& gd, TrajectoryStateOnSurface& tsos) const { + VLocalValues vlp; + ReturnType tuple = getParameters(cluster, gd, tsos); + vlp.emplace_back(std::get<0>(tuple), std::get<1>(tuple)); + return vlp; + } + + + MTDClusterParameterEstimator() {}; + +}; + +#endif diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterizerBase.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterizerBase.h new file mode 100644 index 0000000000000..c665f07aaa46c --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDClusterizerBase.h @@ -0,0 +1,73 @@ +#ifndef RecoLocalTracker_MTDClusterizer_MTDClusterizerBase_H +#define RecoLocalTracker_MTDClusterizer_MTDClusterizerBase_H + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" + +#include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h" +#include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h" + +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" + +#include +#include + + +/** + * Abstract interface for MTD Clusterizers + */ +class MTDClusterizerBase { +public: + + typedef FTLRecHitCollection::const_iterator RecHitIterator; + typedef FTLClusterCollection::const_iterator ClusterIterator; + + // Virtual destructor, this is a base class. + virtual ~MTDClusterizerBase() {} + + // Build clusters + virtual void clusterize( const FTLRecHitCollection & input, + const MTDGeometry* geom, + const MTDTopology* topo, + FTLClusterCollection& output) = 0; + + protected: + struct AccretionCluster { + typedef unsigned short UShort; + static constexpr UShort MAXSIZE = 256; + + std::array energy; + std::array time; + std::array timeError; + std::array x; + std::array y; + + UShort xmin=16000; + UShort ymin=16000; + unsigned int isize=0; + unsigned int curr=0; + + // stack interface (unsafe ok for use below) + UShort top() const { return curr;} + void pop() { ++curr;} + bool empty() { return curr==isize;} + + bool add(FTLCluster::FTLHitPos const & p, float const ienergy, float const itime, float const itimeError) { + if (isize==MAXSIZE) return false; + xmin=std::min(xmin,(unsigned short)(p.row())); + ymin=std::min(ymin,(unsigned short)(p.col())); + energy[isize]=ienergy; + time[isize]=itime; + timeError[isize]=itimeError; + x[isize]=p.row(); + y[isize]=p.col(); + isize++; + return true; + } + }; + + +}; + +#endif diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h new file mode 100644 index 0000000000000..390355f9d3aa4 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h @@ -0,0 +1,91 @@ +#ifndef RecoLocalTracker_FTLClusterizer_MTDThresholdClusterizer_H +#define RecoLocalTracker_FTLClusterizer_MTDThresholdClusterizer_H + +//----------------------------------------------------------------------- +//! \class MTDThresholdClusterizer +//! \brief An explicit threshold-based clustering algorithm. +//! +//! A threshold-based clustering algorithm which clusters FTLRecHits +//! into FTLClusters for each DetUnit. The algorithm is straightforward +//! and purely topological: the clustering process starts with seed hits +//! and continues by adding adjacent hits above the hit threshold. +//! Once the cluster is made, it has to be above the cluster threshold as +//! well. +//! +//! The clusterization is performed on a matrix with size +//! equal to the size of the MTD detector, each cell containing +//! the cahrge and time of the corresponding hit +//! The matrix is reset after each clusterization. +//! +//! The search starts from seed hits, i.e. hits with sufficiently +//! large amplitudes +//! +//! FTLCluster contains a barycenter, but it should be noted that that +//! information is largely useless. One must use a PositionEstimator +//! class to compute the RecHit position and its error for every given +//! cluster. +//! +//! Sets the MTDArrayBuffer dimensions and pixel thresholds. +//! Makes clusters and stores them in theCache if the option +//! useCache has been set. +//----------------------------------------------------------------------- + +// Base class, defines FTLRecHit and FLTCluster. The latter includes +// FTLHit, FTLHitPos and Shift as inner classes. +// +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterizerBase.h" + +// The private pixel buffer +#include "MTDArrayBuffer.h" + +// Parameter Set: +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include + +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" + +class MTDThresholdClusterizer : public MTDClusterizerBase { + public: + + MTDThresholdClusterizer(edm::ParameterSet const& conf); + ~MTDThresholdClusterizer() override; + + // Full I/O in DetSet + void clusterize( const FTLRecHitCollection & input, + const MTDGeometry * geom, + const MTDTopology * topo, + FTLClusterCollection& output) override; + + static void fillDescriptions(edm::ParameterSetDescription& desc); + + private: + + std::vector theSeeds; // cached seed pixels + std::vector theClusters; // resulting clusters + + //! Clustering-related quantities: + float theHitThreshold; // Hit threshold + float theSeedThreshold; // MTD cluster seed + float theClusterThreshold; // Cluster threshold + + //! Geometry-related information + int theNumOfRows; + int theNumOfCols; + + DetId theCurrentId; + + //! Data storage + MTDArrayBuffer theBuffer; // internal nrow * ncol matrix + bool bufferAlreadySet; // status of the buffer array + + bool setup(const MTDGeometry * geometry, const MTDTopology * topo, const DetId& id); + void copy_to_buffer( RecHitIterator itr); + void clear_buffer( RecHitIterator itr); + FTLCluster make_cluster( const FTLCluster::FTLHitPos& hit ); +}; + +#endif diff --git a/RecoLocalFastTime/FTLClusterizer/plugins/BuildFile.xml b/RecoLocalFastTime/FTLClusterizer/plugins/BuildFile.xml new file mode 100644 index 0000000000000..980e858b53f65 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/plugins/BuildFile.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/RecoLocalFastTime/FTLClusterizer/plugins/MTDCPEESProducer.cc b/RecoLocalFastTime/FTLClusterizer/plugins/MTDCPEESProducer.cc new file mode 100644 index 0000000000000..c98cd2c731b27 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/plugins/MTDCPEESProducer.cc @@ -0,0 +1,63 @@ +#include "RecoLocalFastTime/Records/interface/MTDCPERecord.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h" + +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" + +#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/ParameterSet/interface/ParameterSet.h" + +#include +#include + +using namespace edm; + +class MTDCPEESProducer: public edm::ESProducer +{ + public: + MTDCPEESProducer(const edm::ParameterSet & p); + ~MTDCPEESProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); + + std::unique_ptr produce(const MTDCPERecord &); + + private: + edm::ParameterSet pset_; +}; + + +MTDCPEESProducer::MTDCPEESProducer(const edm::ParameterSet & p) +{ + pset_ = p; + setWhatProduced(this,"MTDCPEBase"); +} + +// Configuration descriptions +void +MTDCPEESProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + descriptions.add("MTDCPEESProducer", desc); +} + +MTDCPEESProducer::~MTDCPEESProducer() {} + +std::unique_ptr +MTDCPEESProducer::produce(const MTDCPERecord & iRecord) +{ + edm::ESHandle pDD; + iRecord.getRecord().get( pDD ); + + return std::make_unique( + pset_, + *pDD.product() + ); +} + +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +DEFINE_FWK_EVENTSETUP_MODULE(MTDCPEESProducer); diff --git a/RecoLocalFastTime/FTLClusterizer/plugins/MTDClusterProducer.cc b/RecoLocalFastTime/FTLClusterizer/plugins/MTDClusterProducer.cc new file mode 100644 index 0000000000000..5376054a56755 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/plugins/MTDClusterProducer.cc @@ -0,0 +1,178 @@ +//--------------------------------------------------------------------------- +//! \class MTDClusterProducer +//! +//! \brief EDProducer to cluster FTLRecHits into FTLClusters. +//! +//--------------------------------------------------------------------------- +// Our own stuff +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterizerBase.h" + +// Data Formats +#include "DataFormats/FTLRecHit/interface/FTLRecHit.h" + +// STL +#include +#include +#include +#include + +// MessageLogger +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h" +#include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h" +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" +#include "Geometry/Records/interface/MTDTopologyRcd.h" +#include "Geometry/CommonTopologies/interface/Topology.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/Utilities/interface/InputTag.h" + +class MTDClusterProducer : public edm::stream::EDProducer<> { + public: + //--- Constructor, virtual destructor (just in case) + explicit MTDClusterProducer(const edm::ParameterSet& conf); + ~MTDClusterProducer() override; + static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); + + //--- The top-level event method. + void produce(edm::Event& e, const edm::EventSetup& c) override; + + + //--- Execute the algorithm(s). + template + void run(const T& input, + FTLClusterCollection & output); + + private: + const edm::EDGetTokenT< FTLRecHitCollection > btlHits_; + const edm::EDGetTokenT< FTLRecHitCollection > etlHits_; + + const std::string ftlbInstance_; // instance name of barrel clusters + const std::string ftleInstance_; // instance name of endcap clusters + + const std::string clusterMode_; // user's choice of the clusterizer + std::unique_ptr clusterizer_; // what we got (for now, one ptr to base class) + + const MTDGeometry* geom_; + const MTDTopology* topo_; +}; + + +//--------------------------------------------------------------------------- +//! Constructor: set the ParameterSet and defer all thinking to setupClusterizer(). +//--------------------------------------------------------------------------- +MTDClusterProducer::MTDClusterProducer(edm::ParameterSet const& conf) + : + btlHits_(consumes< FTLRecHitCollection >( conf.getParameter("srcBarrel"))), + etlHits_(consumes< FTLRecHitCollection >( conf.getParameter("srcEndcap"))), + ftlbInstance_(conf.getParameter("BarrelClusterName")), + ftleInstance_(conf.getParameter("EndcapClusterName")), + clusterMode_( conf.getParameter("ClusterMode") ), + clusterizer_(nullptr) // the default, in case we fail to make one +{ + + //--- Declare to the EDM what kind of collections we will be making. + produces(ftlbInstance_); + produces(ftleInstance_); + + //--- Make the algorithm(s) according to what the user specified + //--- in the ParameterSet. + if ( clusterMode_ == "MTDThresholdClusterizer" ) { + clusterizer_ = std::make_unique(conf); + } + else { + throw cms::Exception("MTDClusterProducer") << "[MTDClusterProducer]:" + <<" choice " << clusterMode_ << " is invalid.\n" + << "Possible choices:\n" + << " MTDThresholdClusterizer"; + } +} + +// Destructor +MTDClusterProducer::~MTDClusterProducer() { +} + +// Configuration descriptions +void +MTDClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("srcBarrel", edm::InputTag("mtdRecHits:FTLBarrel")); + desc.add("srcEndcap", edm::InputTag("mtdRecHits:FTLEndcap")); + desc.add("BarrelClusterName", "FTLBarrel"); + desc.add("EndcapClusterName", "FTLEndcap"); + desc.add("ClusterMode","MTDThresholdClusterizer"); + desc.add("HitThreshold", 0.); + desc.add("SeedThreshold", 0.); + desc.add("ClusterThreshold", 0.); + descriptions.add("mtdClusterProducer", desc); +} + +//--------------------------------------------------------------------------- +//! The "Event" entrypoint: gets called by framework for every event +//--------------------------------------------------------------------------- +void MTDClusterProducer::produce(edm::Event& e, const edm::EventSetup& es) +{ + // Step A.1: get input data + edm::Handle< FTLRecHitCollection > inputBarrel; + edm::Handle< FTLRecHitCollection > inputEndcap; + e.getByToken(btlHits_, inputBarrel); + e.getByToken(etlHits_, inputEndcap); + + // Step A.2: get event setup + edm::ESHandle geom; + es.get().get(geom); + geom_ = geom.product(); + + edm::ESHandle mtdTopo; + es.get().get(mtdTopo); + topo_ = mtdTopo.product(); + + // Step B: create the final output collection + auto outputBarrel = std::make_unique(); + auto outputEndcap = std::make_unique(); + + run(*inputBarrel, *outputBarrel ); + run(*inputEndcap, *outputEndcap ); + + e.put(std::move(outputBarrel),ftlbInstance_); + e.put(std::move(outputEndcap),ftleInstance_); +} + + + +//--------------------------------------------------------------------------- +//! Iterate over DetUnits, and invoke the PixelClusterizer on each. +//--------------------------------------------------------------------------- +template +void MTDClusterProducer::run(const T& input, + FTLClusterCollection & output) +{ + if ( ! clusterizer_ ) { + throw cms::Exception("MTDClusterProducer") + <<" at least one clusterizer is not ready -- can't run!" ; + } + + clusterizer_->clusterize( input , geom_, topo_, output); + + LogDebug ("MTDClusterProducer") << " Executing " + << clusterMode_ << " resulted in " << output.size() + << " MTDClusters for " << input.size() << " Hits."; +} + + + + +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +DEFINE_FWK_MODULE(MTDClusterProducer); + diff --git a/RecoLocalFastTime/FTLClusterizer/python/MTDCPEESProducers_cff.py b/RecoLocalFastTime/FTLClusterizer/python/MTDCPEESProducers_cff.py new file mode 100644 index 0000000000000..d9eb9a941c2df --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/python/MTDCPEESProducers_cff.py @@ -0,0 +1,3 @@ +import FWCore.ParameterSet.Config as cms + +from RecoLocalFastTime.FTLClusterizer.MTDCPEESProducer_cfi import * diff --git a/RecoLocalFastTime/FTLClusterizer/python/mtdClusters_cfi.py b/RecoLocalFastTime/FTLClusterizer/python/mtdClusters_cfi.py new file mode 100644 index 0000000000000..f5a5b8e6ac362 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/python/mtdClusters_cfi.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms +from RecoLocalFastTime.FTLClusterizer.mtdClusterProducer_cfi import mtdClusterProducer + +mtdClusters = mtdClusterProducer.clone() diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc new file mode 100644 index 0000000000000..53b6e93a6b500 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc @@ -0,0 +1,114 @@ +#include "Geometry/MTDGeometryBuilder/interface/MTDGeomDetUnit.h" +#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" +#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" + +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h" + +// MessageLogger +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +// Magnetic field +#include "MagneticField/Engine/interface/MagneticField.h" + + +#include + +using namespace std; + +//----------------------------------------------------------------------------- +// A constructor run for generic and templates +// +//----------------------------------------------------------------------------- +MTDCPEBase::MTDCPEBase(edm::ParameterSet const & conf, + const MTDGeometry& geom) + : geom_(geom) +{ + fillDetParams(); +} + + +//----------------------------------------------------------------------------- +// Fill all variables which are constant for an event (geometry) +//----------------------------------------------------------------------------- +void MTDCPEBase::fillDetParams() +{ + auto const & dus = geom_.detUnits(); + unsigned detectors = dus.size(); + m_DetParams.resize(detectors); + LogDebug("MTDCPEBase::fillDetParams():") <<"caching "<(dus[i]); + assert(p.theDet); + + p.theOrigin = p.theDet->surface().toLocal(GlobalPoint(0,0,0)); + + //--- p.theDet->type() returns a GeomDetType, which implements subDetector() + p.thePart = p.theDet->type().subDetector(); + + //--- bounds() is implemented in BoundSurface itself. + p.theThickness = p.theDet->surface().bounds().thickness(); + + // Cache the det id for templates and generic erros + p.theTopol = &(static_cast(p.theDet->topology())); + assert(p.theTopol); + p.theRecTopol = &(static_cast(p.theTopol->specificTopology())); + assert(p.theRecTopol); + + //--- The geometrical description of one module/plaquette + std::pair pitchxy = p.theRecTopol->pitch(); + p.thePitchX = pitchxy.first; // pitch along x + p.thePitchY = pitchxy.second; // pitch along y + + LogDebug("MTDCPEBase::fillDetParams()") << "***** MTD LAYOUT *****" + << " thePart = " << p.thePart + << " theThickness = " << p.theThickness + << " thePitchX = " << p.thePitchX + << " thePitchY = " << p.thePitchY; + } +} + +//----------------------------------------------------------------------------- +// One function to cache the variables common for one DetUnit. +//----------------------------------------------------------------------------- +void +MTDCPEBase::setTheClu( DetParam const & dp, ClusterParam & cp ) const +{ +} + +//------------------------------------------------------------------------ +MTDCPEBase::DetParam const & MTDCPEBase::detParam(const GeomDetUnit & det) const +{ + return m_DetParams.at(det.index()); +} + +LocalPoint +MTDCPEBase::localPosition(DetParam const & dp, ClusterParam & cp) const +{ + //remember measurement point is row(col)+0.5f + MeasurementPoint pos(cp.theCluster->x(),cp.theCluster->y()); + return dp.theTopol->localPosition(pos); +} + +LocalError +MTDCPEBase::localError(DetParam const & dp, ClusterParam & cp) const +{ + constexpr double one_over_twelve = 1./12.; + MeasurementPoint pos(cp.theCluster->x(),cp.theCluster->y()); + MeasurementError simpleRect(one_over_twelve,0,one_over_twelve); + return dp.theTopol->localError(pos,simpleRect); +} + +MTDCPEBase::TimeValue +MTDCPEBase::clusterTime(DetParam const & dp, ClusterParam & cp) const +{ + return cp.theCluster->time(); +} + + +MTDCPEBase::TimeValueError +MTDCPEBase::clusterTimeError(DetParam const & dp, ClusterParam & cp) const +{ + return cp.theCluster->timeError(); +} diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc new file mode 100644 index 0000000000000..c8dc3c358eebe --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -0,0 +1,276 @@ +// Our own includes +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h" + +#include "DataFormats/ForwardDetId/interface/BTLDetId.h" +#include "DataFormats/ForwardDetId/interface/ETLDetId.h" + +#include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h" +#include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h" + +// STL +#include +#include +#include +#include +using namespace std; + +//#define DEBUG_ENABLED +#ifdef DEBUG_ENABLED +#define DEBUG(x) do { std::cout << x << std::endl; } while (0) +#else +#define DEBUG(x) +#endif + +//---------------------------------------------------------------------------- +//! Constructor: +//---------------------------------------------------------------------------- +MTDThresholdClusterizer::MTDThresholdClusterizer + (edm::ParameterSet const& conf) : + // Get energy thresholds + theHitThreshold( conf.getParameter("HitThreshold") ), + theSeedThreshold( conf.getParameter("SeedThreshold") ), + theClusterThreshold( conf.getParameter("ClusterThreshold") ), + theNumOfRows(0), theNumOfCols(0), theCurrentId(0), + theBuffer(theNumOfRows, theNumOfCols ), + bufferAlreadySet(false) +{ +} +///////////////////////////////////////////////////////////////////////////// +MTDThresholdClusterizer::~MTDThresholdClusterizer() {} + + +// Configuration descriptions +void +MTDThresholdClusterizer::fillDescriptions(edm::ParameterSetDescription& desc) { + desc.add("HitThreshold", 0.); + desc.add("SeedThreshold", 0.); + desc.add("ClusterThreshold", 0.); +} + +//---------------------------------------------------------------------------- +//! Prepare the Clusterizer to work on a particular DetUnit. Re-init the +//! size of the panel/plaquette (so update nrows and ncols), +//---------------------------------------------------------------------------- +bool MTDThresholdClusterizer::setup(const MTDGeometry* geom, const MTDTopology* topo, const DetId& id) +{ + theCurrentId=id; + //using geopraphicalId here + const auto& thedet = geom->idToDet(id); + if( thedet == nullptr ) { + throw cms::Exception("MTDThresholdClusterizer") << "GeographicalID: " << std::hex + << id.rawId() + << " is invalid!" << std::dec + << std::endl; + return false; + } + const ProxyMTDTopology& topoproxy = static_cast(thedet->topology()); + const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); + + // Get the new sizes. + unsigned int nrows = topol.nrows(); // rows in x + unsigned int ncols = topol.ncolumns(); // cols in y + + theNumOfRows = nrows; // Set new sizes + theNumOfCols = ncols; + + DEBUG("Buffer size [" << theNumOfRows << "," << theNumOfCols << "]"); + + if ( nrows > theBuffer.rows() || + ncols > theBuffer.columns() ) + { // change only when a larger is needed + // Resize the buffer + theBuffer.setSize(nrows,ncols); // Modify + bufferAlreadySet = true; + } + + return true; +} +//---------------------------------------------------------------------------- +//! \brief Cluster hits. +//! This method operates on a matrix of hits +//! and finds the largest contiguous cluster around +//! each seed hit. +//! Input and output data stored in DetSet +//---------------------------------------------------------------------------- +void MTDThresholdClusterizer::clusterize( const FTLRecHitCollection & input, + const MTDGeometry* geom, + const MTDTopology* topo, + FTLClusterCollection& output) { + + FTLRecHitCollection::const_iterator begin = input.begin(); + FTLRecHitCollection::const_iterator end = input.end(); + + // Do not bother for empty detectors + if (begin == end) + { + edm::LogInfo("MTDThresholdClusterizer") << "No hits to clusterize"; + return; + } + + DEBUG("Input collection " << input.size()); + assert(output.empty()); + + std::set geoIds; + std::multimap geoIdToIdx; + + unsigned index = 0; + for(const auto& hit : input) + { + MTDDetId mtdId=MTDDetId(hit.detid()); + if (mtdId.subDetector() != MTDDetId::FastTime) + { + throw cms::Exception("MTDThresholdClusterizer") << "MTDDetId: " << std::hex + << mtdId.rawId() + << " is invalid!" << std::dec + << std::endl; + } + + if ( mtdId.mtdSubDetector() == MTDDetId::BTL ) + { + BTLDetId hitId(hit.detid()); + DetId geoId = hitId.geographicalId( (BTLDetId::CrysLayout) topo->getMTDTopologyMode() ); //for BTL topology gives different layout id + geoIdToIdx.emplace(geoId,index); + geoIds.emplace(geoId); + ++index; + } + else if ( mtdId.mtdSubDetector() == MTDDetId::ETL ) + { + ETLDetId hitId(hit.detid()); + DetId geoId = hitId.geographicalId(); + geoIdToIdx.emplace(geoId,index); + geoIds.emplace(geoId); + ++index; + } + else + { + throw cms::Exception("MTDThresholdClusterizer") << "MTDDetId: " << std::hex + << mtdId.rawId() + << " is invalid!" << std::dec + << std::endl; + } + } + + //cluster hits within geoIds (modules) + for(unsigned id : geoIds) { + // Set up the clusterization on this DetId. + if ( !setup(geom,topo,DetId(id)) ) + return; + + auto range = geoIdToIdx.equal_range(id); + DEBUG("Matching Ids for " << std::hex << id << std::dec << " [" << range.first->second << "," << range.second->second << "]"); + + // Copy MTDRecHits to the buffer array; select the seed hits + // on the way, and store them in theSeeds. + for(auto itr = range.first; itr != range.second; ++itr) { + const unsigned hitidx = itr->second; + copy_to_buffer(begin+hitidx); + } + + FTLClusterCollection::FastFiller clustersOnDet(output,id); + + for (unsigned int i = 0; i < theSeeds.size(); i++) + { + if ( theBuffer.energy(theSeeds[i]) > theSeedThreshold ) + { // Is this seed still valid? + // Make a cluster around this seed + const FTLCluster & cluster = make_cluster( theSeeds[i] ); + + // Check if the cluster is above threshold + if ( cluster.energy() > theClusterThreshold) + { + DEBUG("putting in this cluster " << i << " #hits:" << cluster.size() + << " E:" << cluster.energy() + << " T:" << cluster.time() + << " X:" << cluster.x() + << " Y:" << cluster.y()); + clustersOnDet.push_back( cluster ); + } + } + } + + // Erase the seeds. + theSeeds.clear(); + // Need to clean unused hits from the buffer array. + for(auto itr = range.first; itr != range.second; ++itr) { + const unsigned hitidx = itr->second; + clear_buffer(begin+hitidx); + } + } +} + +//---------------------------------------------------------------------------- +//! \brief Clear the internal buffer array. +//! +//! MTDs which are not part of recognized clusters are NOT ERASED +//! during the cluster finding. Erase them now. +//! +//---------------------------------------------------------------------------- +void MTDThresholdClusterizer::clear_buffer( RecHitIterator itr ) +{ + theBuffer.clear( itr->row(), itr->column() ); +} + +//---------------------------------------------------------------------------- +//! \brief Copy FTLRecHit into the buffer, identify seeds. +//---------------------------------------------------------------------------- +void MTDThresholdClusterizer::copy_to_buffer( RecHitIterator itr ) +{ + int row = itr->row(); + int col = itr->column(); + float energy = itr->energy(); + float time = itr->time(); + float timeError = itr->timeError(); + + DEBUG("ROW " << row << " COL " << col << " ENERGY " << energy << " TIME " << time); + if ( energy > theHitThreshold) { + theBuffer.set( row, col, energy , time, timeError); + if ( energy > theSeedThreshold) theSeeds.push_back( FTLCluster::FTLHitPos(row,col)); + //sort seeds? + } +} + + +//---------------------------------------------------------------------------- +//! \brief The actual clustering algorithm: group the neighboring hits around the seed. +//---------------------------------------------------------------------------- +FTLCluster +MTDThresholdClusterizer::make_cluster( const FTLCluster::FTLHitPos& hit ) +{ + + //First we acquire the seeds for the clusters + float seed_energy= theBuffer.energy(hit.row(), hit.col()); + float seed_time= theBuffer.time(hit.row(), hit.col()); + float seed_time_error= theBuffer.time_error(hit.row(), hit.col()); + theBuffer.clear(hit); + + AccretionCluster acluster; + acluster.add(hit, seed_energy, seed_time, seed_time_error); + + bool stopClus=false; + //Here we search all hits adjacent to all hits in the cluster. + while ( ! acluster.empty() && ! stopClus) + { + //This is the standard algorithm to find and add a hit + auto curInd = acluster.top(); acluster.pop(); + for ( auto c = std::max(0,int(acluster.y[curInd]-1)); c < std::min(int(acluster.y[curInd]+2),int(theBuffer.columns())) && !stopClus; ++c) { + for ( auto r = std::max(0,int(acluster.x[curInd]-1)); r < std::min(int(acluster.x[curInd]+2),int(theBuffer.rows())) && !stopClus; ++r) { + if ( theBuffer.energy(r,c) > theHitThreshold) { + FTLCluster::FTLHitPos newhit(r,c); + if (!acluster.add( newhit, theBuffer.energy(r,c), theBuffer.time(r,c), theBuffer.time_error(r,c))) + { + stopClus=true; + break; + } + theBuffer.clear(newhit); + } + } + } + } // while accretion + + FTLCluster cluster( theCurrentId, acluster.isize, + acluster.energy.data(), acluster.time.data(), acluster.timeError.data(), + acluster.x.data(), acluster.y.data(), + acluster.xmin, acluster.ymin); + return cluster; +} diff --git a/RecoLocalFastTime/FTLClusterizer/src/module.cc b/RecoLocalFastTime/FTLClusterizer/src/module.cc new file mode 100644 index 0000000000000..e899af8d475d5 --- /dev/null +++ b/RecoLocalFastTime/FTLClusterizer/src/module.cc @@ -0,0 +1,8 @@ +#include "FWCore/Utilities/interface/typelookup.h" + +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h" + + +//--- Now use the Framework macros to set it all up: +// +TYPELOOKUP_DATA_REG(MTDClusterParameterEstimator); diff --git a/RecoLocalFastTime/FTLRecProducers/BuildFile.xml b/RecoLocalFastTime/FTLRecProducers/BuildFile.xml index 6d3f1cab7d6bc..e67e485914b6b 100644 --- a/RecoLocalFastTime/FTLRecProducers/BuildFile.xml +++ b/RecoLocalFastTime/FTLRecProducers/BuildFile.xml @@ -1,7 +1,7 @@ - + diff --git a/RecoLocalFastTime/FTLRecProducers/plugins/BuildFile.xml b/RecoLocalFastTime/FTLRecProducers/plugins/BuildFile.xml index 491b45bb84679..f627694d01de8 100644 --- a/RecoLocalFastTime/FTLRecProducers/plugins/BuildFile.xml +++ b/RecoLocalFastTime/FTLRecProducers/plugins/BuildFile.xml @@ -1,4 +1,6 @@ + + diff --git a/RecoLocalFastTime/FTLRecProducers/plugins/MTDRecHitProducer.cc b/RecoLocalFastTime/FTLRecProducers/plugins/MTDRecHitProducer.cc index ea530ffad3078..7e7dd15f37ab5 100644 --- a/RecoLocalFastTime/FTLRecProducers/plugins/MTDRecHitProducer.cc +++ b/RecoLocalFastTime/FTLRecProducers/plugins/MTDRecHitProducer.cc @@ -9,6 +9,14 @@ #include "RecoLocalFastTime/FTLCommonAlgos/interface/MTDRecHitAlgoBase.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" +#include "Geometry/CommonTopologies/interface/Topology.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" + #include "FWCore/MessageLogger/interface/MessageLogger.h" class MTDRecHitProducer : public edm::stream::EDProducer<> { @@ -27,6 +35,9 @@ class MTDRecHitProducer : public edm::stream::EDProducer<> { const std::string ftleInstance_; // instance name of endcap hits std::unique_ptr barrel_,endcap_; + + edm::ESWatcher geomwatcher_; + const MTDGeometry* geom_; }; MTDRecHitProducer::MTDRecHitProducer(const edm::ParameterSet& ps) : @@ -57,6 +68,9 @@ MTDRecHitProducer::~MTDRecHitProducer() { void MTDRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) { + edm::ESHandle geom; + es.get().get(geom); + geom_ = geom.product(); // tranparently get things from event setup barrel_->getEventSetup(es); @@ -74,7 +88,7 @@ MTDRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) { barrelRechits->reserve(hBarrel->size()/2); for(const auto& uhit : *hBarrel) { uint32_t flags = FTLRecHit::kGood; - auto rechit = std::move( barrel_->makeRecHit(uhit, flags) ); + auto rechit = barrel_->makeRecHit(uhit, flags); if( flags == FTLRecHit::kGood ) barrelRechits->push_back( std::move(rechit) ); } @@ -83,13 +97,14 @@ MTDRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) { endcapRechits->reserve(hEndcap->size()/2); for(const auto& uhit : *hEndcap) { uint32_t flags = FTLRecHit::kGood; - auto rechit = std::move( endcap_->makeRecHit(uhit, flags) ); + auto rechit = endcap_->makeRecHit(uhit, flags); if( flags == FTLRecHit::kGood ) endcapRechits->push_back( std::move(rechit) ); - } - + } + // put the collection of recunstructed hits in the event - evt.put(std::move(barrelRechits), ftlbInstance_); - evt.put(std::move(endcapRechits), ftleInstance_); + // get the orphan handles so we can make refs for the tracking rechits + auto barrelHandle = evt.put(std::move(barrelRechits), ftlbInstance_); + auto endcapHandle = evt.put(std::move(endcapRechits), ftleInstance_); } #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/RecoLocalFastTime/FTLRecProducers/plugins/MTDTrackingRecHitProducer.cc b/RecoLocalFastTime/FTLRecProducers/plugins/MTDTrackingRecHitProducer.cc new file mode 100644 index 0000000000000..a0b5794597fc4 --- /dev/null +++ b/RecoLocalFastTime/FTLRecProducers/plugins/MTDTrackingRecHitProducer.cc @@ -0,0 +1,154 @@ +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h" + +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" +#include "Geometry/CommonTopologies/interface/Topology.h" +#include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h" +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" + +#include "RecoLocalFastTime/Records/interface/MTDCPERecord.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDClusterParameterEstimator.h" +#include "RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +//#define DEBUG_ENABLED +#ifdef DEBUG_ENABLED +#define DEBUG(x) do { std::cout << x << std::endl; } while (0) +#else +#define DEBUG(x) +#endif + +class MTDTrackingRecHitProducer : public edm::stream::EDProducer<> { + + public: + explicit MTDTrackingRecHitProducer(const edm::ParameterSet& ps); + ~MTDTrackingRecHitProducer() override; + static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); + + void produce(edm::Event& evt, const edm::EventSetup& es) override; + void run(edm::Handle > inputHandle, + MTDTrackingDetSetVector &output); + + private: + const edm::EDGetTokenT ftlbClusters_; // collection of barrel digis + const edm::EDGetTokenT ftleClusters_; // collection of endcap digis + + const MTDGeometry* geom_; + const MTDClusterParameterEstimator* cpe_; +}; + +MTDTrackingRecHitProducer::MTDTrackingRecHitProducer(const edm::ParameterSet& ps) : + ftlbClusters_( consumes( ps.getParameter("barrelClusters") ) ), + ftleClusters_( consumes( ps.getParameter("endcapClusters") ) ) +{ + produces< MTDTrackingDetSetVector >(); +} + +// Configuration descriptions +void +MTDTrackingRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("barrelClusters", edm::InputTag("mtdClusters:FTLBarrel")); + desc.add("endcapClusters", edm::InputTag("mtdClusters:FTLEndcap")); + descriptions.add("mtdTrackingRecHitProducer", desc); +} + +MTDTrackingRecHitProducer::~MTDTrackingRecHitProducer() { +} + +void +MTDTrackingRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) { + + edm::ESHandle geom; + es.get().get(geom); + geom_ = geom.product(); + + edm::ESHandle cpe; + es.get().get("MTDCPEBase",cpe); + cpe_ = cpe.product(); + + if ( ! geom_ ) + { + throw cms::Exception("MTDTrackingRecHitProducer") << "Geometry is not available -- can't run!"; + return; // clusterizer is invalid, bail out + } + + if ( ! cpe_ ) + { + throw cms::Exception("MTDTrackingRecHitProducer") << "CPE is not ready -- can't run!"; + return; // clusterizer is invalid, bail out + } + + edm::Handle< FTLClusterCollection > inputBarrel; + evt.getByToken( ftlbClusters_, inputBarrel); + + edm::Handle< FTLClusterCollection > inputEndcap; + evt.getByToken( ftleClusters_, inputEndcap); + + auto outputhits = std::make_unique(); + auto& theoutputhits = *outputhits; + + run(inputBarrel,theoutputhits); + run(inputEndcap,theoutputhits); + + evt.put(std::move(outputhits)); +} + +//--------------------------------------------------------------------------- +//! Iterate over DetUnits, then over Clusters and invoke the CPE on each, +//! and make a RecHit to store the result. +//--------------------------------------------------------------------------- +void MTDTrackingRecHitProducer::run(edm::Handle inputHandle, + MTDTrackingDetSetVector &output) +{ + const edmNew::DetSetVector& input = *inputHandle; + edmNew::DetSetVector::const_iterator DSViter=input.begin(); + + DEBUG("inputCollection " << input.size()); + for ( ; DSViter != input.end() ; DSViter++) { + unsigned int detid = DSViter->detId(); + DetId detIdObject( detid ); + const auto& genericDet = geom_->idToDetUnit(detIdObject); + if( genericDet == nullptr ) { + throw cms::Exception("MTDTrackingRecHitProducer") << "GeographicalID: " << std::hex + << detid + << " is invalid!" << std::dec + << std::endl; + } + + MTDTrackingDetSetVector::FastFiller recHitsOnDet(output,detid); + + edmNew::DetSet::const_iterator clustIt = DSViter->begin(), clustEnd = DSViter->end(); + + for ( ; clustIt != clustEnd; clustIt++) { + DEBUG("Cluster: size " << clustIt->size() << " " << clustIt->x() << "," << clustIt->y() << " " << clustIt->energy() << " " << clustIt->time()); + MTDClusterParameterEstimator::ReturnType tuple = cpe_->getParameters( *clustIt, *genericDet ); + LocalPoint lp( std::get<0>(tuple) ); + LocalError le( std::get<1>(tuple) ); + + + // Create a persistent edm::Ref to the cluster + edm::Ref< edmNew::DetSetVector, FTLCluster > cluster = edmNew::makeRefTo( inputHandle, clustIt); + // Make a RecHit and add it to the DetSet + MTDTrackingRecHit hit( lp, le, *genericDet, cluster); + DEBUG("MTD_TRH: " << hit.localPosition().x() << "," << hit.localPosition().y() + << " : " << hit.localPositionError().xx() << "," << hit.localPositionError().yy() + << " : " << hit.time() << " : " << hit.timeError()); + // Now save it ================= + recHitsOnDet.push_back(hit); + } // <-- End loop on Clusters + } // <-- End loop on DetUnits + DEBUG("outputCollection " << output.size()); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE( MTDTrackingRecHitProducer ); diff --git a/RecoLocalFastTime/FTLRecProducers/python/mtdTrackingRecHits_cfi.py b/RecoLocalFastTime/FTLRecProducers/python/mtdTrackingRecHits_cfi.py new file mode 100644 index 0000000000000..a0f1368adcf32 --- /dev/null +++ b/RecoLocalFastTime/FTLRecProducers/python/mtdTrackingRecHits_cfi.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms +from RecoLocalFastTime.FTLRecProducers.mtdTrackingRecHitProducer_cfi import mtdTrackingRecHitProducer + +mtdTrackingRecHits = mtdTrackingRecHitProducer.clone() diff --git a/RecoLocalFastTime/Records/BuildFile.xml b/RecoLocalFastTime/Records/BuildFile.xml new file mode 100644 index 0000000000000..8730c2cad2bb1 --- /dev/null +++ b/RecoLocalFastTime/Records/BuildFile.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/RecoLocalFastTime/Records/interface/MTDCPERecord.h b/RecoLocalFastTime/Records/interface/MTDCPERecord.h new file mode 100644 index 0000000000000..b25501584c949 --- /dev/null +++ b/RecoLocalFastTime/Records/interface/MTDCPERecord.h @@ -0,0 +1,13 @@ +#ifndef RecoLocalFastTime_Records_MTDCPERecord_h +#define RecoLocalFastTime_Records_MTDCPERecord_h + +#include "FWCore/Framework/interface/EventSetupRecordImplementation.h" +#include "FWCore/Framework/interface/DependentRecordImplementation.h" +#include "Geometry/Records/interface/MTDDigiGeometryRecord.h" + +#include "boost/mpl/vector.hpp" + +class MTDCPERecord: public edm::eventsetup::DependentRecordImplementation > {}; + +#endif diff --git a/RecoLocalFastTime/Records/src/MTDCPERecord.cc b/RecoLocalFastTime/Records/src/MTDCPERecord.cc new file mode 100644 index 0000000000000..f48a9f56319ff --- /dev/null +++ b/RecoLocalFastTime/Records/src/MTDCPERecord.cc @@ -0,0 +1,4 @@ +#include "RecoLocalFastTime/Records/interface/MTDCPERecord.h" +#include "FWCore/Framework/interface/eventsetuprecord_registration_macro.h" + +EVENTSETUP_RECORD_REG(MTDCPERecord); diff --git a/RecoMTD/Configuration/python/RecoMTD_EventContent_cff.py b/RecoMTD/Configuration/python/RecoMTD_EventContent_cff.py new file mode 100644 index 0000000000000..4504e6e3d5fca --- /dev/null +++ b/RecoMTD/Configuration/python/RecoMTD_EventContent_cff.py @@ -0,0 +1,12 @@ +import FWCore.ParameterSet.Config as cms + +#FEVT +RecoMTDFEVT = cms.PSet( + outputCommands = cms.untracked.vstring( + 'keep *_trackExtenderWithMTD_*_*' + ) +) +#RECO content +RecoMTDRECO = RecoMTDFEVT.copy() +#AOD content +RecoMTDAOD = RecoMTDFEVT.copy() diff --git a/RecoMTD/Configuration/python/RecoMTD_cff.py b/RecoMTD/Configuration/python/RecoMTD_cff.py new file mode 100644 index 0000000000000..f1d790f644ecb --- /dev/null +++ b/RecoMTD/Configuration/python/RecoMTD_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +from RecoMTD.TrackExtender.trackExtenderWithMTD_cfi import * + +fastTimingGlobalReco = cms.Sequence(trackExtenderWithMTD) diff --git a/RecoMTD/DetLayers/src/MTDDetLayerGeometry.cc b/RecoMTD/DetLayers/src/MTDDetLayerGeometry.cc index 5d94d78710009..80bec213ae2ff 100644 --- a/RecoMTD/DetLayers/src/MTDDetLayerGeometry.cc +++ b/RecoMTD/DetLayers/src/MTDDetLayerGeometry.cc @@ -105,7 +105,7 @@ MTDDetLayerGeometry::allLayers() const { const DetLayer* MTDDetLayerGeometry::idToLayer(const DetId &id) const{ DetId idout; - MTDDetId detId; + MTDDetId detId(id); if(detId.mtdSubDetector() == 2){ // 2 is ETL ETLDetId etlId( detId.rawId() ); diff --git a/RecoMTD/MeasurementDet/BuildFile.xml b/RecoMTD/MeasurementDet/BuildFile.xml new file mode 100644 index 0000000000000..65c30e468541d --- /dev/null +++ b/RecoMTD/MeasurementDet/BuildFile.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/RecoMTD/MeasurementDet/interface/MTDDetLayerMeasurements.h b/RecoMTD/MeasurementDet/interface/MTDDetLayerMeasurements.h new file mode 100644 index 0000000000000..76f952c222a90 --- /dev/null +++ b/RecoMTD/MeasurementDet/interface/MTDDetLayerMeasurements.h @@ -0,0 +1,128 @@ +#ifndef MeasurementDet_MTDDetLayerMeasurements_H +#define MeasurementDet_MTDDetLayerMeasurements_H + +/** \class MTDDetLayerMeasurements + * The class to access recHits and TrajectoryMeasurements from DetLayer. + * + * \author B. Tannenwald + * Adapted from RecoMuon version. + * + */ + +#include "FWCore/Framework/interface/Event.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" +#include "TrackingTools/TransientTrackingRecHit/interface/GenericTransientTrackingRecHit.h" +#include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h" +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" + +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" + + + +#include + +class DetLayer; +class GeomDet; +class TrajectoryMeasurement; + +class MTDDetLayerMeasurements { + public: + typedef std::vector MeasurementContainer; + typedef std::pair DetWithState; + typedef std::vector MTDRecHitContainer; + + MTDDetLayerMeasurements(edm::InputTag mtdlabel, + edm::ConsumesCollector& iC); + + virtual ~MTDDetLayerMeasurements(); + + // for a given det and state. Not clear when the fastMeasurements below + // should be used, since it isn't passed a GeomDet + MeasurementContainer + measurements( const DetLayer* layer, + const GeomDet * det, + const TrajectoryStateOnSurface& stateOnDet, + const MeasurementEstimator& est, + const edm::Event& iEvent); + + /// returns TMeasurements in a DetLayer compatible with the TSOS. + MeasurementContainer + measurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + const edm::Event& iEvent); + + /// faster version in case the TrajectoryState on the surface of the GeomDet is already available + MeasurementContainer + fastMeasurements( const DetLayer* layer, + const TrajectoryStateOnSurface& theStateOnDet, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + const edm::Event& iEvent); + + /// returns TMeasurements in a DetLayer compatible with the TSOS. + MeasurementContainer + measurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est); + + /// faster version in case the TrajectoryState on the surface of the GeomDet is already available + MeasurementContainer + fastMeasurements( const DetLayer* layer, + const TrajectoryStateOnSurface& theStateOnDet, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est); + + std::vector + groupedMeasurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + const edm::Event& iEvent); + + std::vector + groupedMeasurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est); + + void setEvent(const edm::Event &); + + /// returns the rechits which are on the layer + MTDRecHitContainer recHits(const DetLayer* layer, const edm::Event& iEvent); + + /// returns the rechits which are on the layer + MTDRecHitContainer recHits(const DetLayer* layer); + + + private: + + /// obtain TrackingRecHits from a DetLayer + MTDRecHitContainer recHits(const GeomDet*, const edm::Event& iEvent); + + /// check that the event is set, and throw otherwise + void checkEvent() const; + + + edm::EDGetTokenT mtdToken_; + + // caches that should get filled once per event + edm::Handle> theMTDRecHits; + + void checkMTDRecHits(); + + // keeps track of which event the cache holds + edm::Event::CacheIdentifier_t theMTDEventCacheID; + + const edm::Event* theEvent; +}; +#endif + diff --git a/RecoMTD/MeasurementDet/src/MTDDetLayerMeasurements.cc b/RecoMTD/MeasurementDet/src/MTDDetLayerMeasurements.cc new file mode 100644 index 0000000000000..07f6595f45b38 --- /dev/null +++ b/RecoMTD/MeasurementDet/src/MTDDetLayerMeasurements.cc @@ -0,0 +1,249 @@ +/** \class MTDDetLayerMeasurements + * The class to access recHits and TrajectoryMeasurements from DetLayer. + * + * \author B. Tannenwald + * Adapted from RecoMuon version. + * + */ + +#include "RecoMTD/MeasurementDet/interface/MTDDetLayerMeasurements.h" + +#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" + +#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +typedef std::shared_ptr MTDRecHitPointer; +typedef std::vector MTDRecHitContainer; +typedef MTDDetLayerMeasurements::MeasurementContainer MeasurementContainer; + + +MTDDetLayerMeasurements::MTDDetLayerMeasurements(edm::InputTag mtdlabel, + edm::ConsumesCollector& iC): + theMTDRecHits(), + theMTDEventCacheID(0), + theEvent(nullptr) +{ + mtdToken_ = iC.consumes(mtdlabel); +} + +MTDDetLayerMeasurements::~MTDDetLayerMeasurements(){} + +MTDRecHitContainer MTDDetLayerMeasurements::recHits(const GeomDet* geomDet, + const edm::Event& iEvent) +{ + DetId geoId = geomDet->geographicalId(); + theEvent = &iEvent; + MTDRecHitContainer result; + + checkMTDRecHits(); + + // Create the ChamberId + DetId detId(geoId.rawId()); + LogDebug("Track|RecoMTD|MTDDetLayerMeasurements") << "(MTD): "< bool { return one < two; }; + auto detset = (*theMTDRecHits)[detId]; + + for (const auto & rechit : detset) + result.push_back(GenericTransientTrackingRecHit::build(geomDet, &rechit)); + + return result; +} + +void MTDDetLayerMeasurements::checkMTDRecHits() +{ + LogDebug("Track|RecoMTD|MTDDetLayerMeasurements") << "Checking MTD RecHits"; + checkEvent(); + auto cacheID = theEvent->cacheIdentifier(); + if (cacheID == theMTDEventCacheID) return; + + { + theEvent->getByToken(mtdToken_, theMTDRecHits); + theMTDEventCacheID = cacheID; + } + if(!theMTDRecHits.isValid()) + { + throw cms::Exception("MTDDetLayerMeasurements") << "Cannot get MTD RecHits"; + LogDebug("Track|RecoMTD|MTDDetLayerMeasurements") << "Cannot get MTD RecHits"; + } +} + +///measurements method if already got the Event +MeasurementContainer +MTDDetLayerMeasurements::measurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) { + checkEvent(); + return measurements(layer, startingState, prop, est, *theEvent); +} + + +MeasurementContainer +MTDDetLayerMeasurements::measurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + const edm::Event& iEvent) { + + MeasurementContainer result; + + std::vector dss = layer->compatibleDets(startingState, prop, est); + LogDebug("RecoMTD")<<"compatibleDets: "<dimension() + <<" Chi2: "< +MTDDetLayerMeasurements::groupedMeasurements( const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est) { + checkEvent(); + return groupedMeasurements(layer, startingState, prop, est, *theEvent); +} + + +std::vector +MTDDetLayerMeasurements::groupedMeasurements(const DetLayer* layer, + const TrajectoryStateOnSurface& startingState, + const Propagator& prop, + const MeasurementEstimator& est, + const edm::Event& iEvent) { + + std::vector result; + // if we want to use the concept of InvalidRecHits, + // we can reuse LayerMeasurements from TrackingTools/MeasurementDet + std::vector groups(layer->groupedCompatibleDets(startingState, prop, est)); + + for (const auto& grp : groups) { + + std::vector groupMeasurements; + for (const auto& detAndStateItr : grp) { + + std::vector detMeasurements + = measurements(layer, detAndStateItr.det(), detAndStateItr.trajectoryState(), est, iEvent); + groupMeasurements.insert(groupMeasurements.end(), detMeasurements.begin(), detMeasurements.end()); + } + + if (!groupMeasurements.empty()) + std::sort( groupMeasurements.begin(), groupMeasurements.end(), TrajMeasLessEstim()); + + result.push_back(TrajectoryMeasurementGroup(groupMeasurements, grp)); + } + + return result; +} + +///set event +void MTDDetLayerMeasurements::setEvent(const edm::Event& event) { + theEvent = &event; +} + + +void MTDDetLayerMeasurements::checkEvent() const { + if(!theEvent) + throw cms::Exception("MTDDetLayerMeasurements") << "The event has not been set"; +} + +MTDRecHitContainer MTDDetLayerMeasurements::recHits(const DetLayer* layer, + const edm::Event& iEvent) { + MTDRecHitContainer rhs; + + std::vector gds = layer->basicComponents(); + + for (const GeomDet* igd : gds) { + MTDRecHitContainer detHits = recHits(igd, iEvent); + rhs.insert(rhs.end(), detHits.begin(), detHits.end()); + } + return rhs; +} + +MTDRecHitContainer MTDDetLayerMeasurements::recHits(const DetLayer* layer) +{ + checkEvent(); + return recHits(layer, *theEvent); +} + diff --git a/RecoMTD/TrackExtender/BuildFile.xml b/RecoMTD/TrackExtender/BuildFile.xml new file mode 100644 index 0000000000000..76cfe56b18cd2 --- /dev/null +++ b/RecoMTD/TrackExtender/BuildFile.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/RecoMTD/TrackExtender/plugins/BuildFile.xml b/RecoMTD/TrackExtender/plugins/BuildFile.xml new file mode 100644 index 0000000000000..99fc4330d8d6b --- /dev/null +++ b/RecoMTD/TrackExtender/plugins/BuildFile.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc b/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc new file mode 100644 index 0000000000000..4465d28fb9987 --- /dev/null +++ b/RecoMTD/TrackExtender/plugins/TrackExtenderWithMTD.cc @@ -0,0 +1,610 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h" +#include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" + +#include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" + +#include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h" +#include "RecoMTD/DetLayers/interface/MTDDetTray.h" +#include "RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h" +#include "RecoMTD/DetLayers/interface/MTDDetRing.h" + +#include "DataFormats/ForwardDetId/interface/BTLDetId.h" +#include "DataFormats/ForwardDetId/interface/ETLDetId.h" +#include "DataFormats/ForwardDetId/interface/MTDChannelIdentifier.h" +#include "Geometry/CommonTopologies/interface/PixelTopology.h" + +#include "TrackingTools/PatternTools/interface/Trajectory.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" + +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" + +#include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHitBuilder.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" + +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" + +#include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h" + +#include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h" +#include "TrackingTools/TrackRefitter/interface/TrackTransformer.h" + +#include + +#include "Geometry/CommonTopologies/interface/Topology.h" + +using namespace std; +using namespace edm; + + + + +template +class TrackExtenderWithMTDT : public edm::stream::EDProducer<> { + static constexpr char pathLengthName[] = "pathLength"; + static constexpr char pathLengthOrigTrkName[] = "generalTrackPathLength"; + static constexpr char betaOrigTrkName[] = "generalTrackBeta"; + static constexpr char t0OrigTrkName[] = "generalTrackt0"; + static constexpr char covt0t0OrigTrkName[] = "generalTrackcovt0t0"; + public: + typedef typename TrackCollection:: value_type TrackType; + typedef edm::View InputCollection; + + TrackExtenderWithMTDT(const ParameterSet& pset); + + void produce(edm::Event& ev, const edm::EventSetup& es) final; + + TransientTrackingRecHit::ConstRecHitContainer tryBTLLayers(const TrackType&, + const MTDTrackingDetSetVector&, + const MTDDetLayerGeometry*, + const MagneticField* field, + const Propagator* prop) const; + + TransientTrackingRecHit::ConstRecHitContainer tryETLLayers(const TrackType&, + const MTDTrackingDetSetVector&, + const MTDDetLayerGeometry*, + const MagneticField* field, + const Propagator* prop) const; + + RefitDirection::GeometricalDirection + checkRecHitsOrdering(TransientTrackingRecHit::ConstRecHitContainer const & recHits) const { + + if (!recHits.empty()){ + GlobalPoint first = gtg->idToDet(recHits.front()->geographicalId())->position(); + GlobalPoint last = gtg->idToDet(recHits.back()->geographicalId())->position(); + + // maybe perp2? + auto rFirst = first.mag2(); + auto rLast = last.mag2(); + if(rFirst < rLast) return RefitDirection::insideOut; + if(rFirst > rLast) return RefitDirection::outsideIn; + } + LogDebug("TrackExtenderWithMTD") << "Impossible to determine the rechits order" < tracksToken_; + edm::EDGetTokenT hitsToken_; + edm::EDGetTokenT bsToken_; + const bool updateTraj_, updateExtra_, updatePattern_; + const std::string mtdRecHitBuilder_,propagator_, transientTrackBuilder_; + std::unique_ptr theEstimator; + std::unique_ptr theTransformer; + edm::ESHandle builder; + edm::ESHandle hitbuilder; + edm::ESHandle gtg; + edm::ESHandle prop; +}; + + +template +TrackExtenderWithMTDT::TrackExtenderWithMTDT(const ParameterSet& iConfig) : + tracksToken_(consumes(iConfig.getParameter("tracksSrc"))), + hitsToken_(consumes(iConfig.getParameter("hitsSrc"))), + bsToken_(consumes(iConfig.getParameter("beamSpotSrc"))), + updateTraj_(iConfig.getParameter("updateTrackTrajectory")), + updateExtra_(iConfig.getParameter("updateTrackExtra")), + updatePattern_(iConfig.getParameter("updateTrackHitPattern")), + mtdRecHitBuilder_(iConfig.getParameter("MTDRecHitBuilder")), + propagator_(iConfig.getParameter("Propagator")), + transientTrackBuilder_(iConfig.getParameter("TransientTrackBuilder")) { + constexpr float maxChi2=25.; + constexpr float nSigma=3.; + theEstimator = std::make_unique(maxChi2,nSigma); + + theTransformer = std::make_unique(iConfig.getParameterSet("TrackTransformer")); + + produces >(pathLengthOrigTrkName); + produces >(betaOrigTrkName); + produces >(t0OrigTrkName); + produces >(covt0t0OrigTrkName); + produces >(pathLengthName); + produces>(); + produces(); + produces(); +} + +template +void TrackExtenderWithMTDT::produce( edm::Event& ev, + const edm::EventSetup& es ) { + //this produces pieces of the track extra + Traj2TrackHits t2t; + + theTransformer->setServices(es); + + TrackingRecHitRefProd hitsRefProd = ev.getRefBeforePut(); + reco::TrackExtraRefProd extrasRefProd = ev.getRefBeforePut(); + + es.get().get(gtg); + + edm::ESHandle geo; + es.get().get(geo); + + edm::ESHandle magfield; + es.get().get(magfield); + + es.get().get(transientTrackBuilder_, builder); + es.get().get(mtdRecHitBuilder_,hitbuilder); + + edm::ESHandle prop; + es.get().get(propagator_,prop); + + edm::ESHandle httopo; + es.get().get(httopo); + const TrackerTopology& ttopo = *httopo; + + auto output = std::make_unique(); + auto extras = std::make_unique(); + auto outhits = std::make_unique>(); + + auto pathLengths = std::make_unique>(); + std::vector pathLengthsRaw; + + auto pathLengthsOrigTrk = std::make_unique>(); + std::vector pathLengthsOrigTrkRaw; + + auto betaOrigTrk = std::make_unique>(); + std::vector betaOrigTrkRaw; + + auto t0OrigTrk = std::make_unique>(); + std::vector t0OrigTrkRaw; + + auto covt0t0OrigTrk = std::make_unique>(); + std::vector covt0t0OrigTrkRaw; + + edm::Handle tracksH; + ev.getByToken(tracksToken_,tracksH); + const auto& tracks = *tracksH; + + edm::Handle hitsH; + ev.getByToken(hitsToken_,hitsH); + const auto& hits = *hitsH; + + edm::Handle bsH; + ev.getByToken(bsToken_,bsH); + const auto& bs = *bsH; + + std::vector track_indices; + unsigned itrack = 0; + for( const auto& track : tracks ) { + + reco::TransientTrack ttrack(track,magfield.product(),gtg); + auto trajs = theTransformer->transform(track); + auto thits = theTransformer->getTransientRecHits(ttrack); + + TransientTrackingRecHit::ConstRecHitContainer mtdthits; + for( auto& ahit : tryBTLLayers(track,hits,geo.product(),magfield.product(),prop.product()) ) { + mtdthits.push_back(ahit); + } + // in the future this should include an intermediate refit before propagating to the ETL + // for now it is ok + for( auto& ahit : tryETLLayers(track,hits,geo.product(),magfield.product(),prop.product()) ) { + mtdthits.push_back(ahit); + } + + auto ordering = checkRecHitsOrdering(thits); + if( ordering == RefitDirection::insideOut) { + for( auto& ahit : mtdthits ) thits.push_back(ahit); + } else { + std::reverse(mtdthits.begin(),mtdthits.end()); + for( auto& ahit : thits ) mtdthits.push_back(ahit); + thits.swap(mtdthits); + } + auto trajwithmtd = theTransformer->transform(ttrack,thits); + float pathLengthMap = -1.f, betaMap = 0.f, t0Map = 0.f, covt0t0Map = -1.f; + + for( const auto& trj : trajwithmtd ) { + const auto& thetrj = (updateTraj_ ? trj : trajs.front()); + float pathLength = 0.f; + reco::Track result = buildTrack(track, thetrj, trj, bs, magfield.product(), + prop.product(), !mtdthits.empty(),pathLength); + if( result.ndof() >= 0 ) { + /// setup the track extras + reco::TrackExtra::TrajParams trajParams; + reco::TrackExtra::Chi2sFive chi2s; + size_t hitsstart = outhits->size(); + if( updatePattern_ ) { + t2t(trj,*outhits,trajParams,chi2s); // this fills the output hit collection + } else { + t2t(thetrj,*outhits,trajParams,chi2s); + } + size_t hitsend = outhits->size(); + extras->push_back(buildTrackExtra(trj)); // always push back the fully built extra, update by setting in track + extras->back().setHits(hitsRefProd,hitsstart,hitsend-hitsstart); + extras->back().setTrajParams(trajParams,chi2s); + //create the track + output->push_back(result); + pathLengthsRaw.push_back(pathLength); + pathLengthMap = pathLength; + auto& backtrack = output->back(); + betaMap = backtrack.beta(); + t0Map = backtrack.t0(); + covt0t0Map = backtrack.covt0t0(); + reco::TrackExtraRef extraRef(extrasRefProd,extras->size()-1); + backtrack.setExtra( (updateExtra_ ? extraRef : track.extra()) ); + for(unsigned ihit = hitsstart; ihit < hitsend; ++ihit) { + backtrack.appendHitPattern((*outhits)[ihit],ttopo); + } + } + } + pathLengthsOrigTrkRaw.push_back(pathLengthMap); + betaOrigTrkRaw.push_back(betaMap); + t0OrigTrkRaw.push_back(t0Map); + covt0t0OrigTrkRaw.push_back(covt0t0Map); + ++itrack; + } + + auto outTrksHandle = ev.put(std::move(output)); + ev.put(std::move(extras)); + ev.put(std::move(outhits)); + + edm::ValueMap::Filler fillerPathLengths(*pathLengths); + fillerPathLengths.insert(outTrksHandle,pathLengthsRaw.cbegin(),pathLengthsRaw.cend()); + fillerPathLengths.fill(); + ev.put(std::move(pathLengths),pathLengthName); + + edm::ValueMap::Filler fillerPathLengthsOrigTrk(*pathLengthsOrigTrk); + fillerPathLengthsOrigTrk.insert(tracksH,pathLengthsOrigTrkRaw.cbegin(),pathLengthsOrigTrkRaw.cend()); + fillerPathLengthsOrigTrk.fill(); + ev.put(std::move(pathLengthsOrigTrk),pathLengthOrigTrkName); + + edm::ValueMap::Filler fillerBetas(*betaOrigTrk); + fillerBetas.insert(tracksH,betaOrigTrkRaw.cbegin(),betaOrigTrkRaw.cend()); + fillerBetas.fill(); + ev.put(std::move(betaOrigTrk),betaOrigTrkName); + + edm::ValueMap::Filler fillert0s(*t0OrigTrk); + fillert0s.insert(tracksH,t0OrigTrkRaw.cbegin(),t0OrigTrkRaw.cend()); + fillert0s.fill(); + ev.put(std::move(t0OrigTrk),t0OrigTrkName); + + edm::ValueMap::Filler fillercovt0t0s(*covt0t0OrigTrk); + fillercovt0t0s.insert(tracksH,covt0t0OrigTrkRaw.cbegin(),covt0t0OrigTrkRaw.cend()); + fillercovt0t0s.fill(); + ev.put(std::move(covt0t0OrigTrk),covt0t0OrigTrkName); +} + +namespace { + auto cmp = [](const unsigned one, const unsigned two) -> bool { return one < two; }; + + void find_hits_in_dets(const MTDTrackingDetSetVector& hits, const DetLayer* layer, + const TrajectoryStateOnSurface& tsos, const Propagator* prop, + const MeasurementEstimator& theEstimator, + const TransientTrackingRecHitBuilder& hitbuilder, + TransientTrackingRecHit::ConstRecHitContainer& output) { + pair comp = layer->compatible(tsos,*prop,theEstimator); + if( comp.first ) { + vector compDets = layer->compatibleDets(tsos,*prop,theEstimator); + if (!compDets.empty()) { + for( const auto& detWithState : compDets ) { + auto range = hits.equal_range(detWithState.first->geographicalId(),cmp); + for( auto detitr = range.first; detitr != range.second; ++detitr ) { + auto best = detitr->end(); + double best_chi2 = std::numeric_limits::max(); + for( auto itr = detitr->begin(); itr != detitr->end(); ++itr ) { + auto est = theEstimator.estimate(detWithState.second,*itr); + if( est.first && est.second < best_chi2 ) { // just take the best chi2 + best = itr; + best_chi2 = est.second; + } + } + if( best != detitr->end() ) { + output.push_back(hitbuilder.build(&*best)); + } + } + } + } + } + } +} + +template +TransientTrackingRecHit::ConstRecHitContainer +TrackExtenderWithMTDT::tryBTLLayers(const TrackType& track, + const MTDTrackingDetSetVector& hits, + const MTDDetLayerGeometry* geo, + const MagneticField* field, + const Propagator* prop) const { + std::unique_ptr trkProp(prop->clone()); + + TransientTrackingRecHit::ConstRecHitContainer output; + const vector& layers = geo->allBTLLayers(); + auto tTrack = builder->build(track); + + for (const DetLayer* ilay : layers) { + // get the outermost trajectory point on the track + TrajectoryStateOnSurface tsos = tTrack.outermostMeasurementState(); + find_hits_in_dets(hits,ilay,tsos,trkProp.get(),*theEstimator,*hitbuilder,output); + } + return output; +} + +template +TransientTrackingRecHit::ConstRecHitContainer +TrackExtenderWithMTDT::tryETLLayers(const TrackType& track, + const MTDTrackingDetSetVector& hits, + const MTDDetLayerGeometry* geo, + const MagneticField* field, + const Propagator* prop) const { + std::unique_ptr trkProp(prop->clone()); + + TransientTrackingRecHit::ConstRecHitContainer output; + const vector& layers = geo->allETLLayers(); + + auto tTrack = builder->build(track); + + for (const DetLayer* ilay : layers) { + const BoundDisk& disk = static_cast(ilay)->specificSurface(); + const double diskZ = disk.position().z(); + + // get the outermost trajectory point on the track + TrajectoryStateOnSurface tsos = tTrack.outermostMeasurementState(); + if( tsos.globalPosition().z() * diskZ < 0 ) continue; // only propagate to the disk that's on the same side + find_hits_in_dets(hits,ilay,tsos,trkProp.get(),*theEstimator,*hitbuilder,output); + } + return output; +} + + +//below is unfortunately ripped from other places but +//since track producer doesn't know about MTD we have to do this +template +reco::Track TrackExtenderWithMTDT::buildTrack(const reco::Track& orig, + const Trajectory& traj, + const Trajectory& trajWithMtd, + const reco::BeamSpot& bs, + const MagneticField* field, + const Propagator* thePropagator, + bool hasMTD, + float& pathLengthOut) const { + std::unique_ptr trkProp(thePropagator->clone()); + + // get the state closest to the beamline + TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface = + traj.closestMeasurement(GlobalPoint(bs.x0(),bs.y0(),bs.z0())).updatedState(); + + if UNLIKELY(!stateForProjectionToBeamLineOnSurface.isValid()) { + edm::LogError("CannotPropagateToBeamLine")<<"the state on the closest measurement isnot valid. skipping track."; + return reco::Track(); + } + + constexpr double mpi = 0.13957018; + constexpr double c = 2.99792458e1; //[cm/ns] + + const FreeTrajectoryState & stateForProjectionToBeamLine=*stateForProjectionToBeamLineOnSurface.freeState(); + + TSCBLBuilderWithPropagator tscblBuilder(*trkProp); + TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs); + + if UNLIKELY(!tscbl.isValid()) { + return reco::Track(); + } + + GlobalPoint v = tscbl.trackStateAtPCA().position(); + math::XYZPoint pos( v.x(), v.y(), v.z() ); + GlobalVector p = tscbl.trackStateAtPCA().momentum(); + math::XYZVector mom( p.x(), p.y(), p.z() ); + + int ndof = traj.ndof(); + + double t0 = 0.; + double covt0t0 = -1.; + pathLengthOut = -1.f; // if there is no MTD flag the pathlength with -1 + double betaOut = 0.; + + //compute path length for time backpropagation, using first MTD hit for the momentum + if (hasMTD) { + + bool validpropagation = true; + double pathlength = 0.; + double pathlength1 = 0.; + double pathlength2 = 0.; + for (auto it=trajWithMtd.measurements().begin(); it!=trajWithMtd.measurements().end()-1; ++it) { + const auto &propresult = trkProp->propagateWithPath(it->updatedState(), (it+1)->updatedState().surface()); + double layerpathlength = std::abs(propresult.second); + if (layerpathlength==0.) { + validpropagation = false; + } + pathlength1 += layerpathlength; + } + + double thit = 0.; + double thiterror = -1.; + bool validmtd = false; + if (trajWithMtd.direction() == alongMomentum) { + for (auto it=trajWithMtd.measurements().begin(); it!=trajWithMtd.measurements().end(); ++it) { + bool ismtd = it->recHit()->geographicalId().det() == DetId::Forward && ForwardSubdetector(it->recHit()->geographicalId().subdetId()) == FastTime; + if (ismtd) { + const auto &propresult2 = trkProp->propagateWithPath(tscbl.trackStateAtPCA(), trajWithMtd.firstMeasurement().updatedState().surface()); + pathlength2 = propresult2.second; + if (pathlength2 == 0.) { + validpropagation = false; + } + pathlength = pathlength1 + pathlength2; + const MTDTrackingRecHit *mtdhit = static_cast(it->recHit()->hit()); + thit = mtdhit->time(); + thiterror = mtdhit->timeError(); + validmtd = true; + break; + } + } + } + else { + for (auto it=trajWithMtd.measurements().rbegin(); it!=trajWithMtd.measurements().rend(); ++it) { + bool ismtd = it->recHit()->geographicalId().det() == DetId::Forward && ForwardSubdetector(it->recHit()->geographicalId().subdetId()) == FastTime; + if (ismtd) { + const auto &propresult2 = trkProp->propagateWithPath(tscbl.trackStateAtPCA(), trajWithMtd.lastMeasurement().updatedState().surface()); + pathlength2 = propresult2.second; + if (pathlength2 == 0.) { + validpropagation = false; + } + pathlength = pathlength1 + pathlength2; + const MTDTrackingRecHit *mtdhit = static_cast(it->recHit()->hit()); + thit = mtdhit->time(); + thiterror = mtdhit->timeError(); + validmtd = true; + break; + } + } + } + + if (validmtd && validpropagation) { + + double magp = p.mag(); + double gammasq = 1. + magp*magp/mpi/mpi; + double beta = std::sqrt(1.-1./gammasq); + double dt = pathlength/beta/c; + pathLengthOut = pathlength; // set path length if we've got a timing hit + t0 = thit - dt; + covt0t0 = thiterror*thiterror; + beta = betaOut; + } + } + + return reco::Track(traj.chiSquared(), + int(ndof), + pos, mom, tscbl.trackStateAtPCA().charge(), + tscbl.trackStateAtPCA().curvilinearError(), + orig.algo(),reco::TrackBase::undefQuality,t0,betaOut,covt0t0,-1.); +} + +template +reco::TrackExtra TrackExtenderWithMTDT::buildTrackExtra(const Trajectory& trajectory) const { + + static const string metname = "TrackExtenderWithMTD"; + + const Trajectory::RecHitContainer transRecHits = trajectory.recHits(); + + // put the collection of TrackingRecHit in the event + + // sets the outermost and innermost TSOSs + // ToDo: validation for track states with MTD + TrajectoryStateOnSurface outerTSOS; + TrajectoryStateOnSurface innerTSOS; + unsigned int innerId=0, outerId=0; + TrajectoryMeasurement::ConstRecHitPointer outerRecHit; + DetId outerDetId; + + if (trajectory.direction() == alongMomentum) { + LogTrace(metname)<<"alongMomentum"; + outerTSOS = trajectory.lastMeasurement().updatedState(); + innerTSOS = trajectory.firstMeasurement().updatedState(); + outerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId(); + innerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId(); + outerRecHit = trajectory.lastMeasurement().recHit(); + outerDetId = trajectory.lastMeasurement().recHit()->geographicalId(); + } + else if (trajectory.direction() == oppositeToMomentum) { + LogTrace(metname)<<"oppositeToMomentum"; + outerTSOS = trajectory.firstMeasurement().updatedState(); + innerTSOS = trajectory.lastMeasurement().updatedState(); + outerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId(); + innerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId(); + outerRecHit = trajectory.firstMeasurement().recHit(); + outerDetId = trajectory.firstMeasurement().recHit()->geographicalId(); + } + else LogError(metname)<<"Wrong propagation direction!"; + + const GeomDet *outerDet = gtg->idToDet(outerDetId); + GlobalPoint outerTSOSPos = outerTSOS.globalParameters().position(); + bool inside = outerDet->surface().bounds().inside(outerDet->toLocal(outerTSOSPos)); + + + GlobalPoint hitPos = (outerRecHit->isValid()) ? outerRecHit->globalPosition() : outerTSOS.globalParameters().position() ; + + if(!inside) { + LogTrace(metname) + << "The Global Muon outerMostMeasurementState is not compatible with the recHit detector!" + << " Setting outerMost postition to recHit position if recHit isValid: " << outerRecHit->isValid(); + LogTrace(metname)<<"From " << outerTSOSPos << " to " << hitPos; + } + + + //build the TrackExtra + GlobalPoint v = (inside) ? outerTSOSPos : hitPos ; + GlobalVector p = outerTSOS.globalParameters().momentum(); + math::XYZPoint outpos( v.x(), v.y(), v.z() ); + math::XYZVector outmom( p.x(), p.y(), p.z() ); + + v = innerTSOS.globalParameters().position(); + p = innerTSOS.globalParameters().momentum(); + math::XYZPoint inpos( v.x(), v.y(), v.z() ); + math::XYZVector inmom( p.x(), p.y(), p.z() ); + + reco::TrackExtra trackExtra(outpos, outmom, true, inpos, inmom, true, + outerTSOS.curvilinearError(), outerId, + innerTSOS.curvilinearError(), innerId, + trajectory.direction(),trajectory.seedRef()); + + return trackExtra; + +} + +template +string TrackExtenderWithMTDT::dumpLayer(const DetLayer* layer) const { + stringstream output; + + const BoundSurface* sur=nullptr; + const BoundCylinder* bc=nullptr; + const BoundDisk* bd=nullptr; + + sur = &(layer->surface()); + if ( (bc = dynamic_cast(sur)) ) { + output << " Cylinder of radius: " << bc->radius() << endl; + } + else if ( (bd = dynamic_cast(sur)) ) { + output << " Disk at: " << bd->position().z() << endl; + } + return output.str(); +} + +//define this as a plug-in +#include +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h" +typedef TrackExtenderWithMTDT TrackExtenderWithMTD; + +DEFINE_FWK_MODULE(TrackExtenderWithMTD); diff --git a/RecoMTD/TrackExtender/python/PropagatorWithMaterialForMTD_cfi.py b/RecoMTD/TrackExtender/python/PropagatorWithMaterialForMTD_cfi.py new file mode 100644 index 0000000000000..5e8db4f224dc5 --- /dev/null +++ b/RecoMTD/TrackExtender/python/PropagatorWithMaterialForMTD_cfi.py @@ -0,0 +1,13 @@ +import FWCore.ParameterSet.Config as cms + +PropagatorWithMaterialForMTD = cms.ESProducer( + "PropagatorWithMaterialESProducer", + MaxDPhi = cms.double(1.6), #default was 1.6 + ComponentName = cms.string('PropagatorWithMaterialForMTD'), + Mass = cms.double(0.13957018), #default was 0.105 + PropagationDirection = cms.string('anyDirection'), + useRungeKutta = cms.bool(False), + ptMin = cms.double(0.1), + useOldAnalPropLogic = cms.bool(False) +) + diff --git a/RecoMTD/TrackExtender/python/trackExtenderWithMTD_cfi.py b/RecoMTD/TrackExtender/python/trackExtenderWithMTD_cfi.py new file mode 100644 index 0000000000000..d4623180481a2 --- /dev/null +++ b/RecoMTD/TrackExtender/python/trackExtenderWithMTD_cfi.py @@ -0,0 +1,31 @@ +import FWCore.ParameterSet.Config as cms + +from RecoMTD.TrackExtender.PropagatorWithMaterialForMTD_cfi import * + +_mtdRecHitBuilder = cms.string('MTDRecHitBuilder') +_mtdPropagator = cms.string('PropagatorWithMaterialForMTD') + +trackExtenderWithMTD = cms.EDProducer( + 'TrackExtenderWithMTD', + tracksSrc = cms.InputTag("generalTracks"), + hitsSrc = cms.InputTag("mtdTrackingRecHits"), + beamSpotSrc = cms.InputTag("offlineBeamSpot"), + updateTrackTrajectory = cms.bool(True), + updateTrackExtra = cms.bool(True), + updateTrackHitPattern = cms.bool(True), + TransientTrackBuilder = cms.string('TransientTrackBuilder'), + MTDRecHitBuilder = _mtdRecHitBuilder, + Propagator = _mtdPropagator, + TrackTransformer = cms.PSet( + DoPredictionsOnly = cms.bool(False), + Fitter = cms.string('KFFitterForRefitInsideOut'), + #TrackerRecHitBuilder = cms.string('WithTrackAngleAndTemplate'), + TrackerRecHitBuilder = cms.string('WithTrackAngle'), + Smoother = cms.string('KFSmootherForRefitInsideOut'), + MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = _mtdRecHitBuilder, + RefitDirection = cms.string('alongMomentum'), + RefitRPCHits = cms.bool(True), + Propagator = _mtdPropagator + ) + ) diff --git a/RecoMTD/TransientTrackingRecHit/BuildFile.xml b/RecoMTD/TransientTrackingRecHit/BuildFile.xml new file mode 100644 index 0000000000000..b1908e5bc11c0 --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/BuildFile.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHit.h b/RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHit.h new file mode 100644 index 0000000000000..f54d29acf35c7 --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHit.h @@ -0,0 +1,61 @@ +#ifndef RecoMTD_TransientTrackingRecHit_MTDTransientTrackingRecHit_h +#define RecoMTD_TransientTrackingRecHit_MTDTransientTrackingRecHit_h + +/** \class MTDTransientTrackingRecHit + * + * A TransientTrackingRecHit for MTD. + * + * + * \author L. Gray - FNAL + * + */ + + +#include "TrackingTools/TransientTrackingRecHit/interface/GenericTransientTrackingRecHit.h" +#include "DataFormats/TrackingRecHit/interface/RecSegment.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + + +class MTDTransientTrackingRecHit final : public GenericTransientTrackingRecHit{ +public: + using MTDRecHitPointer = std::shared_ptr; + using ConstMTDRecHitPointer = std::shared_ptr; + + typedef std::vector MTDRecHitContainer; + typedef std::vector ConstMTDRecHitContainer; + + ~MTDTransientTrackingRecHit() override{} + + /// if this rec hit is a BTL rec hit + bool isBTL() const; + + /// if this rec hit is a ETL rec hit + bool isETL() const; + + static RecHitPointer build( const GeomDet * geom, const TrackingRecHit* rh) { + return RecHitPointer( new MTDTransientTrackingRecHit(geom, rh)); + } + + static MTDRecHitPointer specificBuild(const GeomDet * geom, const TrackingRecHit* rh) { + LogDebug("MTDTransientTrackingRecHit") << "Getting specificBuild"< trackingGeometry = nullptr); + + ~MTDTransientTrackingRecHitBuilder() override {} ; + + using TransientTrackingRecHitBuilder::build; + /// Call the MTDTransientTrackingRecHit::specificBuild + RecHitPointer build(const TrackingRecHit *p, + edm::ESHandle trackingGeometry) const ; + + RecHitPointer build(const TrackingRecHit * p) const override; + + ConstRecHitContainer build(const trackingRecHit_iterator& start, const trackingRecHit_iterator& stop) const; + + private: + edm::ESHandle theTrackingGeometry; + +}; + +#endif diff --git a/RecoMTD/TransientTrackingRecHit/plugins/BuildFile.xml b/RecoMTD/TransientTrackingRecHit/plugins/BuildFile.xml new file mode 100644 index 0000000000000..18a3607fcd7a4 --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/plugins/BuildFile.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/RecoMTD/TransientTrackingRecHit/plugins/MTDTransientTrackingRecHitBuilderESProducer.cc b/RecoMTD/TransientTrackingRecHit/plugins/MTDTransientTrackingRecHitBuilderESProducer.cc new file mode 100644 index 0000000000000..e42cee9d616d6 --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/plugins/MTDTransientTrackingRecHitBuilderESProducer.cc @@ -0,0 +1,63 @@ +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" + +#include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHitBuilder.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" + +#include + +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" + +#include "FWCore/Framework/interface/ESProducer.h" + +namespace edm {class ParameterSet;} + +class TransientRecHitRecord; + +class MTDTransientTrackingRecHitBuilderESProducer: public edm::ESProducer { +public: + /// Constructor + MTDTransientTrackingRecHitBuilderESProducer(const edm::ParameterSet&); + + /// Destructor + ~MTDTransientTrackingRecHitBuilderESProducer() override; + + // Operations + std::unique_ptr produce(const TransientRecHitRecord&); + +protected: + +private: +}; + +using namespace edm; +using namespace std; + +MTDTransientTrackingRecHitBuilderESProducer::MTDTransientTrackingRecHitBuilderESProducer(const ParameterSet & parameterSet) { + + setWhatProduced(this,parameterSet.getParameter("ComponentName")); +} + +MTDTransientTrackingRecHitBuilderESProducer::~MTDTransientTrackingRecHitBuilderESProducer() {} + + +std::unique_ptr +MTDTransientTrackingRecHitBuilderESProducer::produce(const TransientRecHitRecord& iRecord){ + + + ESHandle trackingGeometry; + iRecord.getRecord().get(trackingGeometry); + + return std::make_unique(trackingGeometry); +} + +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Utilities/interface/typelookup.h" + +DEFINE_FWK_EVENTSETUP_MODULE(MTDTransientTrackingRecHitBuilderESProducer); diff --git a/RecoMTD/TransientTrackingRecHit/python/MTDTransientTrackingRecHitBuilder_cfi.py b/RecoMTD/TransientTrackingRecHit/python/MTDTransientTrackingRecHitBuilder_cfi.py new file mode 100644 index 0000000000000..9b4a43faaf67e --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/python/MTDTransientTrackingRecHitBuilder_cfi.py @@ -0,0 +1,8 @@ +import FWCore.ParameterSet.Config as cms + +MTDTransientTrackingRecHitBuilder = cms.ESProducer("MTDTransientTrackingRecHitBuilderESProducer", + ComponentName = cms.string('MTDRecHitBuilder') +) + + + diff --git a/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHit.cc b/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHit.cc new file mode 100644 index 0000000000000..a15f04e7c641d --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHit.cc @@ -0,0 +1,40 @@ +/** \file + * + */ + +#include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHit.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" + +#include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h" +#include "DataFormats/ForwardDetId/interface/MTDDetId.h" +#include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "FWCore/Utilities/interface/Exception.h" + +#include + +typedef MTDTransientTrackingRecHit::MTDRecHitPointer MTDRecHitPointer; +typedef MTDTransientTrackingRecHit::RecHitContainer MTDRecHitContainer; + + +MTDTransientTrackingRecHit::MTDTransientTrackingRecHit(const GeomDet* geom, const TrackingRecHit* rh) : + GenericTransientTrackingRecHit(*geom,*rh){} + +MTDTransientTrackingRecHit::MTDTransientTrackingRecHit(const MTDTransientTrackingRecHit& other ) : + GenericTransientTrackingRecHit(*other.det(), *(other.hit())) {} + +bool MTDTransientTrackingRecHit::isBTL() const{ + MTDDetId temp(geographicalId()); + return (temp.mtdSubDetector() == MTDDetId::BTL); +} + +bool MTDTransientTrackingRecHit::isETL() const{ + MTDDetId temp(geographicalId()); + return (temp.mtdSubDetector() == MTDDetId::ETL); +} + +void MTDTransientTrackingRecHit::invalidateHit(){ + setType(bad); //trackingRecHit_->setType(bad); // maybe add in later +} + diff --git a/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHitBuilder.cc b/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHitBuilder.cc new file mode 100644 index 0000000000000..829d0ae42ad34 --- /dev/null +++ b/RecoMTD/TransientTrackingRecHit/src/MTDTransientTrackingRecHitBuilder.cc @@ -0,0 +1,51 @@ +/** + * Class: MTDTransientTrackingRecHitBuilder + * + * Description: + * + * + * + * Authors : + * L. Gray FNAL + * + **/ + +#include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHitBuilder.h" +#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" + +MTDTransientTrackingRecHitBuilder::MTDTransientTrackingRecHitBuilder(edm::ESHandle trackingGeometry): + theTrackingGeometry(trackingGeometry) +{} + + +MTDTransientTrackingRecHitBuilder::RecHitPointer +MTDTransientTrackingRecHitBuilder::build (const TrackingRecHit* p, + edm::ESHandle trackingGeometry) const { + + if ( p->geographicalId().det() == DetId::Forward && + p->geographicalId().subdetId() == FastTime) { + return p->cloneSH(); + } + + return RecHitPointer(); + +} + +MTDTransientTrackingRecHitBuilder::RecHitPointer +MTDTransientTrackingRecHitBuilder::build(const TrackingRecHit * p) const { + if(theTrackingGeometry.isValid()) return build(p,theTrackingGeometry); + else + throw cms::Exception("MTD|RecoMTD|MTDTransientTrackingRecHitBuilder") + <<"ERROR! You are trying to build a MTDTransientTrackingRecHit with a non valid GlobalTrackingGeometry"; +} + +MTDTransientTrackingRecHitBuilder::ConstRecHitContainer +MTDTransientTrackingRecHitBuilder::build(const trackingRecHit_iterator& start, const trackingRecHit_iterator& stop) const { + + ConstRecHitContainer result; + for(trackingRecHit_iterator hit = start; hit != stop; ++hit ) + result.push_back(build(&**hit)); + + return result; +} + diff --git a/RecoMuon/GlobalTrackingTools/python/GlobalTrajectoryBuilderCommon_cff.py b/RecoMuon/GlobalTrackingTools/python/GlobalTrajectoryBuilderCommon_cff.py index be49cef0b514b..f3db179246ef4 100644 --- a/RecoMuon/GlobalTrackingTools/python/GlobalTrajectoryBuilderCommon_cff.py +++ b/RecoMuon/GlobalTrackingTools/python/GlobalTrajectoryBuilderCommon_cff.py @@ -19,6 +19,7 @@ TrackerRecHitBuilder = cms.string('WithAngleAndTemplate'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), DoPredictionsOnly = cms.bool(False) diff --git a/RecoMuon/L3MuonProducer/src/L3MuonProducer.cc b/RecoMuon/L3MuonProducer/src/L3MuonProducer.cc index 91687fbf8e8de..6940cf9649e28 100644 --- a/RecoMuon/L3MuonProducer/src/L3MuonProducer.cc +++ b/RecoMuon/L3MuonProducer/src/L3MuonProducer.cc @@ -269,14 +269,16 @@ void L3MuonProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio psd0.add("PCut", 2.5); { edm::ParameterSetDescription psd1; - psd1.add("DoPredictionsOnly", false); - psd1.add("Fitter", "hltESPL3MuKFTrajectoryFitter"); - psd1.add("TrackerRecHitBuilder", "hltESPTTRHBWithTrackAngle"); - psd1.add("Smoother", "hltESPKFTrajectorySmootherForMuonTrackLoader"); - psd1.add("MuonRecHitBuilder", "hltESPMuonTransientTrackingRecHitBuilder"); - psd1.add("RefitDirection", "insideOut"); - psd1.add("RefitRPCHits", true); - psd1.add("Propagator", "hltESPSmartPropagatorAny"); + TrackTransformer::fillPSetDescription(psd1, + false, // do predictions only + "hltESPL3MuKFTrajectoryFitter", // fitter + "hltESPKFTrajectorySmootherForMuonTrackLoader", // smoother + "hltESPSmartPropagatorAny", // propagator + "insideOut", // refit direction + true, // refit rpc hits + "hltESPTTRHBWithTrackAngle", // tracker rechit builder + "hltESPMuonTransientTrackingRecHitBuilder" // muon rechit builder + ); psd0.add("TrackTransformer", psd1); } { diff --git a/RecoMuon/MuonIdentification/python/TrackerKinkFinder_cfi.py b/RecoMuon/MuonIdentification/python/TrackerKinkFinder_cfi.py index 233ca21d4bd2a..3d511d4eb80ae 100644 --- a/RecoMuon/MuonIdentification/python/TrackerKinkFinder_cfi.py +++ b/RecoMuon/MuonIdentification/python/TrackerKinkFinder_cfi.py @@ -12,6 +12,7 @@ TrackerRecHitBuilder = cms.string('WithAngleAndTemplate'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite'), diff --git a/RecoTracker/SpecialSeedGenerators/python/inOutSeedsFromTrackerMuons_cfi.py b/RecoTracker/SpecialSeedGenerators/python/inOutSeedsFromTrackerMuons_cfi.py index 80fa55c211adf..2826792262d0a 100644 --- a/RecoTracker/SpecialSeedGenerators/python/inOutSeedsFromTrackerMuons_cfi.py +++ b/RecoTracker/SpecialSeedGenerators/python/inOutSeedsFromTrackerMuons_cfi.py @@ -16,6 +16,7 @@ TrackerRecHitBuilder = cms.string('WithAngleAndTemplate'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite'), diff --git a/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc b/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc index ed2d07618a7f9..b016eb91a12a3 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc +++ b/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc @@ -73,8 +73,7 @@ void BTLTileDeviceSim::getHitsResponse(const std::vectoremplace(mtd_digitizer::MTDCellId(id,row,col), mtd_digitizer::MTDCellInfo()).first; diff --git a/TrackingTools/RecoGeometry/BuildFile.xml b/TrackingTools/RecoGeometry/BuildFile.xml index 13948da3c720d..ef807b711dff2 100644 --- a/TrackingTools/RecoGeometry/BuildFile.xml +++ b/TrackingTools/RecoGeometry/BuildFile.xml @@ -2,6 +2,8 @@ + + diff --git a/TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h b/TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h index b89be762bdf84..de7744ed9c0cd 100644 --- a/TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h +++ b/TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h @@ -9,6 +9,7 @@ #include "DataFormats/DetId/interface/DetId.h" #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" +#include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h" #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h" #include "TrackingTools/DetLayers/interface/DetLayerGeometry.h" @@ -18,9 +19,14 @@ class DetLayer; class GlobalDetLayerGeometry: public DetLayerGeometry { public: - GlobalDetLayerGeometry(const GeometricSearchTracker* tracker, + GlobalDetLayerGeometry(const GeometricSearchTracker* tracker, const MuonDetLayerGeometry* muon): - tracker_(tracker),muon_(muon){}; + tracker_(tracker),muon_(muon),mtd_(nullptr){}; + + GlobalDetLayerGeometry(const GeometricSearchTracker* tracker, + const MuonDetLayerGeometry* muon, + const MTDDetLayerGeometry* mtd): + tracker_(tracker),muon_(muon),mtd_(mtd){}; ~GlobalDetLayerGeometry() override {} @@ -38,6 +44,7 @@ class GlobalDetLayerGeometry: public DetLayerGeometry { private: const GeometricSearchTracker* tracker_; const MuonDetLayerGeometry* muon_; + const MTDDetLayerGeometry* mtd_; }; diff --git a/TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h b/TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h index 5064ba29871e6..9d1e427f7bcf6 100644 --- a/TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h +++ b/TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h @@ -6,6 +6,7 @@ #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h" #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h" +#include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h" //#include "RecoTracker/Record/interface/NavigationSchoolRecord.h" //#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" @@ -13,7 +14,7 @@ class RecoGeometryRecord : public edm::eventsetup::DependentRecordImplementation > {}; diff --git a/TrackingTools/RecoGeometry/plugins/DetLayerGeometryESProducer.cc b/TrackingTools/RecoGeometry/plugins/DetLayerGeometryESProducer.cc index 28dd56aab9dd4..79bbd5a0e9693 100644 --- a/TrackingTools/RecoGeometry/plugins/DetLayerGeometryESProducer.cc +++ b/TrackingTools/RecoGeometry/plugins/DetLayerGeometryESProducer.cc @@ -21,8 +21,6 @@ DetLayerGeometryESProducer::~DetLayerGeometryESProducer() {} std::unique_ptr DetLayerGeometryESProducer::produce(const RecoGeometryRecord & iRecord){ - - return std::make_unique(); } diff --git a/TrackingTools/RecoGeometry/plugins/GlobalDetLayerGeometryESProducer.cc b/TrackingTools/RecoGeometry/plugins/GlobalDetLayerGeometryESProducer.cc index 4811812a1fc5f..4d425ea50be15 100644 --- a/TrackingTools/RecoGeometry/plugins/GlobalDetLayerGeometryESProducer.cc +++ b/TrackingTools/RecoGeometry/plugins/GlobalDetLayerGeometryESProducer.cc @@ -23,11 +23,26 @@ GlobalDetLayerGeometryESProducer::produce(const RecoGeometryRecord & iRecord){ edm::ESHandle tracker; edm::ESHandle muon; + edm::ESHandle mtd; iRecord.getRecord().get(tracker); iRecord.getRecord().get(muon); + + // get the MTD if it is available + if(auto mtdRecord = iRecord.tryToGetRecord()) { + mtdRecord->get(mtd); + if(!mtd.isValid()) { + LogInfo("GlobalDetLayergGeometryBuilder") << "No MTD geometry is available."; + } + } else { + LogInfo("GlobalDetLayerGeometryBuilder") << "No MTDDigiGeometryRecord is available."; + } + + // if we've got MTD initialize it + if( mtd.isValid() ) return std::make_unique(tracker.product(), muon.product(), mtd.product()); return std::make_unique(tracker.product(), muon.product()); + } diff --git a/TrackingTools/RecoGeometry/src/GlobalDetLayerGeometry.cc b/TrackingTools/RecoGeometry/src/GlobalDetLayerGeometry.cc index 5b1902e199b0d..e6f2f3011f385 100644 --- a/TrackingTools/RecoGeometry/src/GlobalDetLayerGeometry.cc +++ b/TrackingTools/RecoGeometry/src/GlobalDetLayerGeometry.cc @@ -1,13 +1,18 @@ #include "TrackingTools/RecoGeometry/interface/GlobalDetLayerGeometry.h" +#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" #include "FWCore/Utilities/interface/typelookup.h" const DetLayer* GlobalDetLayerGeometry::idToLayer(const DetId& detId) const{ - if(detId.det() ==1) return tracker_->idToLayer(detId); - else if(detId.det() ==2) return muon_->idToLayer(detId); + if(detId.det() == DetId::Tracker) return tracker_->idToLayer(detId); + else if(detId.det() == DetId::Muon) return muon_->idToLayer(detId); + else if(mtd_ != nullptr && + detId.det() == DetId::Forward && + detId.subdetId() == FastTime ) return mtd_->idToLayer(detId); else{ throw cms::Exception("DetLayers") - << "Error: called GlobalDetLayerGeometry::idToLayer() for a detId which is neither Tracker nor Muon " << detId; + << "Error: called GlobalDetLayerGeometry::idToLayer() for a detId which is neither Tracker nor Muon " << (mtd_ == nullptr ? "" : "nor MTD ") + << detId; } } diff --git a/TrackingTools/TrackRefitter/BuildFile.xml b/TrackingTools/TrackRefitter/BuildFile.xml index 6b1ca8b7398ba..b737c28c30c96 100644 --- a/TrackingTools/TrackRefitter/BuildFile.xml +++ b/TrackingTools/TrackRefitter/BuildFile.xml @@ -4,6 +4,7 @@ + diff --git a/TrackingTools/TrackRefitter/interface/TrackTransformer.h b/TrackingTools/TrackRefitter/interface/TrackTransformer.h index b46d1c4a5550d..0bd7f2c2bc39e 100644 --- a/TrackingTools/TrackRefitter/interface/TrackTransformer.h +++ b/TrackingTools/TrackRefitter/interface/TrackTransformer.h @@ -26,7 +26,7 @@ #include "DataFormats/TrackReco/interface/TrackFwd.h" -namespace edm {class ParameterSet; class EventSetup;} +namespace edm {class ParameterSet; class EventSetup; class ParameterSetDescription;} namespace reco {class TransientTrack;} class TrajectoryFitter; @@ -44,6 +44,18 @@ class TrackTransformer final : public TrackTransformerBase{ /// Destructor ~TrackTransformer() override; + + /// fillDescriptions + static void fillPSetDescription(edm::ParameterSetDescription& descriptions, + bool DoPredictionsOnly=false, + const std::string& Fitter="KFFitterForRefitInsideOut", + const std::string& Smoother="KFSmootherForRefitInsideOut", + const std::string& Propagator="SmartPropagatorAnyRK", + const std::string& RefitDirection="alongMomentum", + bool RefitRPCHits=true, + const std::string& TrackerRecHitBuilder="WithTrackAngle", + const std::string& MuonRecHitBuilder="MuonRecHitBuilder", + const std::string& MTDRecHitBuilder="MTDRecHitBuilder"); // Operations @@ -110,6 +122,10 @@ class TrackTransformer final : public TrackTransformerBase{ const std::string theMuonRecHitBuilderName; edm::ESHandle theMuonRecHitBuilder; + + const std::string theMTDRecHitBuilderName; + bool theMtdAvailable; + edm::ESHandle theMTDRecHitBuilder; }; #endif diff --git a/TrackingTools/TrackRefitter/python/TracksToTrajectories_cff.py b/TrackingTools/TrackRefitter/python/TracksToTrajectories_cff.py index 769f4ddfd664a..b1e2165b4b674 100644 --- a/TrackingTools/TrackRefitter/python/TracksToTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/TracksToTrajectories_cff.py @@ -8,6 +8,7 @@ from TrackingTools.KalmanUpdators.KFUpdatorESProducer_cfi import * from TrackingTools.GeomPropagators.SmartPropagator_cff import * +from RecoMTD.TransientTrackingRecHit.MTDTransientTrackingRecHitBuilder_cfi import * from RecoMuon.TransientTrackingRecHit.MuonTransientTrackingRecHitBuilder_cfi import * from RecoTracker.TransientTrackingRecHit.TransientTrackingRecHitBuilder_cfi import * diff --git a/TrackingTools/TrackRefitter/python/cosmicMuonTrajectories_cff.py b/TrackingTools/TrackRefitter/python/cosmicMuonTrajectories_cff.py index 271088aeeeb9e..63e93440077cb 100644 --- a/TrackingTools/TrackRefitter/python/cosmicMuonTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/cosmicMuonTrajectories_cff.py @@ -25,6 +25,7 @@ Tracks = cms.InputTag("cosmicMuons"), TrackTransformer = cms.PSet( TrackerRecHitBuilder = cms.string('WithTrackAngle'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitRPCHits = cms.bool(True) ) ) diff --git a/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectoriesP5_cff.py b/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectoriesP5_cff.py index 6d3846a347978..f1bacc12fd7c8 100644 --- a/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectoriesP5_cff.py +++ b/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectoriesP5_cff.py @@ -22,6 +22,7 @@ TrackerRecHitBuilder = cms.string('WithTrackAngle'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRK') diff --git a/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectories_cff.py b/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectories_cff.py index 51331b945994e..b7defea421c43 100644 --- a/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/ctfWithMaterialTrajectories_cff.py @@ -22,6 +22,7 @@ TrackerRecHitBuilder = cms.string('WithTrackAngle'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite') diff --git a/TrackingTools/TrackRefitter/python/globalCosmicMuonTrajectories_cff.py b/TrackingTools/TrackRefitter/python/globalCosmicMuonTrajectories_cff.py index 2dbd0afb3ace2..d58b45bc0c071 100644 --- a/TrackingTools/TrackRefitter/python/globalCosmicMuonTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/globalCosmicMuonTrajectories_cff.py @@ -18,6 +18,7 @@ Tracks = cms.InputTag("globalCosmicMuons"), TrackTransformer = cms.PSet(TrackerRecHitBuilder = cms.string('WithTrackAngle'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitRPCHits = cms.bool(True), # muon station to be skipped //also kills RPCs in that station SkipStationDT = cms.int32(-999), diff --git a/TrackingTools/TrackRefitter/python/globalMuonTrajectories_cff.py b/TrackingTools/TrackRefitter/python/globalMuonTrajectories_cff.py index 71eb1dbb2269a..27073dbd5b7ba 100644 --- a/TrackingTools/TrackRefitter/python/globalMuonTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/globalMuonTrajectories_cff.py @@ -22,6 +22,7 @@ TrackerRecHitBuilder = cms.string('WithTrackAngle'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite') diff --git a/TrackingTools/TrackRefitter/python/standAloneMuonTrajectories_cff.py b/TrackingTools/TrackRefitter/python/standAloneMuonTrajectories_cff.py index c6d836f9c7394..f88313bf4c898 100644 --- a/TrackingTools/TrackRefitter/python/standAloneMuonTrajectories_cff.py +++ b/TrackingTools/TrackRefitter/python/standAloneMuonTrajectories_cff.py @@ -22,6 +22,7 @@ TrackerRecHitBuilder = cms.string('WithTrackAngle'), Smoother = cms.string('KFSmootherForRefitInsideOut'), MuonRecHitBuilder = cms.string('MuonRecHitBuilder'), + MTDRecHitBuilder = cms.string('MTDRecHitBuilder'), RefitDirection = cms.string('alongMomentum'), RefitRPCHits = cms.bool(True), Propagator = cms.string('SmartPropagatorAnyRKOpposite') diff --git a/TrackingTools/TrackRefitter/src/TrackTransformer.cc b/TrackingTools/TrackRefitter/src/TrackTransformer.cc index 9e7832ed2ade7..c54362fe46ef2 100644 --- a/TrackingTools/TrackRefitter/src/TrackTransformer.cc +++ b/TrackingTools/TrackRefitter/src/TrackTransformer.cc @@ -22,6 +22,8 @@ #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/DetId/interface/DetId.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + using namespace std; using namespace edm; @@ -34,12 +36,33 @@ TrackTransformer::TrackTransformer(const ParameterSet& parameterSet): theSmootherName(parameterSet.getParameter("Smoother")), thePropagatorName(parameterSet.getParameter("Propagator")), theTrackerRecHitBuilderName(parameterSet.getParameter("TrackerRecHitBuilder")), - theMuonRecHitBuilderName(parameterSet.getParameter("MuonRecHitBuilder")) + theMuonRecHitBuilderName(parameterSet.getParameter("MuonRecHitBuilder")), + theMTDRecHitBuilderName(parameterSet.getParameter("MTDRecHitBuilder")) {} /// Destructor TrackTransformer::~TrackTransformer(){} +void TrackTransformer::fillPSetDescription(edm::ParameterSetDescription& desc, + bool DoPredictionsOnly, + const std::string& Fitter, + const std::string& Smoother, + const std::string& Propagator, + const std::string& RefitDirection, + bool RefitRPCHits, + const std::string& TrackerRecHitBuilder, + const std::string& MuonRecHitBuilder, + const std::string& MTDRecHitBuilder) { + desc.add("DoPredictionsOnly",DoPredictionsOnly); + desc.add("Fitter",Fitter); + desc.add("Smoother",Smoother); + desc.add("Propagator",Propagator); + desc.add("RefitDirection",RefitDirection); + desc.add("RefitRPCHits",RefitRPCHits); + desc.add("TrackerRecHitBuilder",TrackerRecHitBuilder); + desc.add("MuonRecHitBuilder",MuonRecHitBuilder); + desc.add("MTDRecHitBuilder",MTDRecHitBuilder); +} void TrackTransformer::setServices(const EventSetup& setup){ @@ -83,6 +106,8 @@ void TrackTransformer::setServices(const EventSetup& setup){ LogTrace(metname) << "TransientRecHitRecord changed!"; setup.get().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder); setup.get().get(theMuonRecHitBuilderName,theMuonRecHitBuilder); + setup.get().get(theMTDRecHitBuilderName,theMTDRecHitBuilder); + theMtdAvailable = theMTDRecHitBuilder.isValid(); hitCloner = static_cast(theTrackerRecHitBuilder.product())->cloner(); } theFitter->setHitCloner(&hitCloner); @@ -113,6 +138,10 @@ TrackTransformer::getTransientRecHits(const reco::TransientTrack& track) const { continue; } result.push_back(theMuonRecHitBuilder->build(&**hit)); + } else if ( (*hit)->geographicalId().det() == DetId::Forward && + (*hit)->geographicalId().subdetId() == FastTime ) { + if ( theMtdAvailable ) result.push_back(theMTDRecHitBuilder->build(&**hit)); + else throw cms::Exception("TrackTransformer") << "MTD hit encountered but MTD not available!"; } } }