diff --git a/L1Trigger/DTTriggerPhase2/interface/DTprimitive.h b/L1Trigger/DTTriggerPhase2/interface/DTprimitive.h index cec71913409f8..512ffca09cc49 100644 --- a/L1Trigger/DTTriggerPhase2/interface/DTprimitive.h +++ b/L1Trigger/DTTriggerPhase2/interface/DTprimitive.h @@ -37,6 +37,11 @@ class DTPrimitive { const int superLayerId() const { return superLayerId_; }; const cmsdt::LATERAL_CASES laterality() const { return laterality_; }; + bool operator==(const DTPrimitive& dtp) { + return (tdcTimeStamp() == dtp.tdcTimeStamp() && channelId() == dtp.channelId() && layerId() == dtp.layerId() && + cameraId() == dtp.cameraId() && cameraId() == dtp.cameraId() && superLayerId() == dtp.superLayerId()); + } + private: int cameraId_; // Chamber ID int superLayerId_; // SL ID diff --git a/L1Trigger/DTTriggerPhase2/interface/LateralityBasicProvider.h b/L1Trigger/DTTriggerPhase2/interface/LateralityBasicProvider.h new file mode 100644 index 0000000000000..6fbb9d82a4c96 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/LateralityBasicProvider.h @@ -0,0 +1,46 @@ +#ifndef L1Trigger_DTTriggerPhase2_LateralityBasicProvider_h +#define L1Trigger_DTTriggerPhase2_LateralityBasicProvider_h + +#include "L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h" + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +struct lat_combination { + short missing_layer; + short cellLayout[cmsdt::NUM_LAYERS]; + lat_vector latcombs; +}; + +// =============================================================================== +// Class declarations +// =============================================================================== + +class LateralityBasicProvider : public LateralityProvider { +public: + // Constructors and destructor + LateralityBasicProvider(const edm::ParameterSet &pset, edm::ConsumesCollector &iC); + ~LateralityBasicProvider() override; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities) override; + + void finish() override; + + // Other public methods + +private: + // Private methods + void analyze(MuonPathPtr &inMPath, std::vector &lateralities); + void fill_lat_combinations(); + // Private attributes + const bool debug_; + std::vector lat_combinations; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/LateralityCoarsedProvider.h b/L1Trigger/DTTriggerPhase2/interface/LateralityCoarsedProvider.h new file mode 100644 index 0000000000000..e7f0ac25800b2 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/LateralityCoarsedProvider.h @@ -0,0 +1,50 @@ +#ifndef L1Trigger_DTTriggerPhase2_LateralityCoarsedProvider_h +#define L1Trigger_DTTriggerPhase2_LateralityCoarsedProvider_h + +#include "L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h" + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +struct lat_coarsed_combination { + short missing_layer; + short cellLayout[cmsdt::NUM_LAYERS]; + short coarsed_times[cmsdt::NUM_LAYERS]; + lat_vector latcombs; +}; + +// =============================================================================== +// Class declarations +// =============================================================================== + +class LateralityCoarsedProvider : public LateralityProvider { +public: + // Constructors and destructor + LateralityCoarsedProvider(const edm::ParameterSet &pset, edm::ConsumesCollector &iC); + ~LateralityCoarsedProvider() override; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities) override; + + void finish() override; + + // Other public methods + +private: + // Private methods + void analyze(MuonPathPtr &inMPath, std::vector &lateralities); + std::vector coarsify_times(MuonPathPtr &inMPath); + void fill_lat_combinations(); + std::vector> convertString(std::string chain); + // Private attributes + const bool debug_; + std::vector lat_combinations; + edm::FileInPath laterality_filename_; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h b/L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h new file mode 100644 index 0000000000000..866fcb08533cb --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h @@ -0,0 +1,58 @@ +#ifndef Phase2L1Trigger_DTTrigger_LateralityProvider_h +#define Phase2L1Trigger_DTTrigger_LateralityProvider_h + +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/FrameworkfwdMostUsed.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "L1Trigger/DTTriggerPhase2/interface/MuonPath.h" +#include "L1Trigger/DTTriggerPhase2/interface/constants.h" + +#include +#include + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +// =============================================================================== +// Class declarations +// =============================================================================== + +using latcomb = std::vector; +using lat_vector = std::vector; + +class LateralityProvider { +public: + // Constructors and destructor + LateralityProvider(const edm::ParameterSet& pset, edm::ConsumesCollector& iC); + virtual ~LateralityProvider(); + + // Main methods + virtual void initialise(const edm::EventSetup& iEventSetup); + virtual void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + MuonPathPtrs& inMpath, + std::vector& lateralities) = 0; + + virtual void finish(); + + // Other public methods + + // Public attributes + lat_vector LAT_VECTOR_NULL = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}; + +private: + // Private methods + + // Private attributes + const bool debug_; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MPCleanHitsFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPCleanHitsFilter.h index d458ccdce1cea..f3713b322adbd 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MPCleanHitsFilter.h +++ b/L1Trigger/DTTriggerPhase2/interface/MPCleanHitsFilter.h @@ -26,7 +26,11 @@ class MPCleanHitsFilter : public MPFilter { const edm::EventSetup& iEventSetup, std::vector& inMPath, std::vector& outMPath) override{}; - + void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + std::vector& inSLMPath, + std::vector& inCorMPath, + std::vector& outMPath) override{}; void run(edm::Event& iEvent, const edm::EventSetup& iEventSetup, MuonPathPtrs& inMPath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MPCorFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPCorFilter.h new file mode 100644 index 0000000000000..01cd5a860cc69 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MPCorFilter.h @@ -0,0 +1,76 @@ +#ifndef Phase2L1Trigger_DTTrigger_MPCorFilter_h +#define Phase2L1Trigger_DTTrigger_MPCorFilter_h + +#include "L1Trigger/DTTriggerPhase2/interface/MPFilter.h" +#include "L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h" + +#include +#include + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +// =============================================================================== +// Class declarations +// =============================================================================== + +struct valid_cor_tp_t { + bool valid; + cmsdt::metaPrimitive mp; + int coarsed_t0; + int coarsed_pos; + int coarsed_slope; + valid_cor_tp_t() : valid(false), mp(cmsdt::metaPrimitive()), coarsed_t0(-1), coarsed_pos(-1), coarsed_slope(-1) {} + valid_cor_tp_t(bool valid, cmsdt::metaPrimitive mp, int coarsed_t0, int coarsed_pos, int coarsed_slope) + : valid(valid), mp(mp), coarsed_t0(coarsed_t0), coarsed_pos(coarsed_pos), coarsed_slope(coarsed_slope) {} +}; + +using valid_cor_tp_arr_t = std::vector; + +class MPCorFilter : public MPFilter { +public: + // Constructors and destructor + MPCorFilter(const edm::ParameterSet &pset); + ~MPCorFilter() override = default; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPath, + std::vector &outMPath) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inSLMPath, + std::vector &inCorMPath, + std::vector &outMPath) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMPath, + MuonPathPtrs &outMPath) override{}; + + void finish() override; + + // Other public methods + + // Public attributes + void printmP(cmsdt::metaPrimitive mP); + +private: + // Private methods + std::vector filter(std::vector SL1mps, + std::vector SL2mps, + std::vector SL3mps, + std::vector Cormps); + std::vector coarsify(cmsdt::metaPrimitive mp, int sl); + bool isDead(cmsdt::metaPrimitive mp, std::vector coarsed, std::map tps_per_bx); + int killTps(cmsdt::metaPrimitive mp, std::vector coarsed, int bx, std::map &tps_per_bx); + int match(cmsdt::metaPrimitive mp, std::vector coarsed, valid_cor_tp_t valid_cor_tp2); + int get_chi2(cmsdt::metaPrimitive mp); + + // Private attributes + const bool debug_; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MPFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPFilter.h index a75a1a74e6474..c7743bf75f304 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MPFilter.h +++ b/L1Trigger/DTTriggerPhase2/interface/MPFilter.h @@ -12,7 +12,7 @@ #include "Geometry/Records/interface/MuonGeometryRecord.h" #include "Geometry/DTGeometry/interface/DTGeometry.h" -#include "Geometry/DTGeometry/interface/DTLayer.h" +#include "DataFormats/MuonDetId/interface/DTLayerId.h" #include #include @@ -37,6 +37,11 @@ class MPFilter { const edm::EventSetup& iEventSetup, std::vector& inMPath, std::vector& outMPath) = 0; + virtual void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + std::vector& inSLMPath, + std::vector& inCorMPath, + std::vector& outMPath) = 0; virtual void run(edm::Event& iEvent, const edm::EventSetup& iEventSetup, MuonPathPtrs& inMPath, @@ -46,6 +51,12 @@ class MPFilter { // Other public methods + // Public attributes + // max drift velocity + edm::FileInPath maxdrift_filename_; + int maxdriftinfo_[5][4][14]; + int max_drift_tdc = -1; + private: // Private attributes const bool debug_; diff --git a/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilter.h index 03f52b181db2d..03db99c9a6d9b 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilter.h +++ b/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilter.h @@ -26,6 +26,11 @@ class MPQualityEnhancerFilter : public MPFilter { const edm::EventSetup &iEventSetup, std::vector &inMPath, std::vector &outMPath) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inSLMPath, + std::vector &inCorMPath, + std::vector &outMPath) override{}; void run(edm::Event &iEvent, const edm::EventSetup &iEventSetup, MuonPathPtrs &inMPath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilterBayes.h b/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilterBayes.h index 5db4357905d5a..ef108561bb06b 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilterBayes.h +++ b/L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilterBayes.h @@ -27,6 +27,11 @@ class MPQualityEnhancerFilterBayes : public MPFilter { const edm::EventSetup &iEventSetup, std::vector &inMPath, std::vector &outMPath) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inSLMPath, + std::vector &inCorMPath, + std::vector &outMPath) override{}; void run(edm::Event &iEvent, const edm::EventSetup &iEventSetup, MuonPathPtrs &inMPath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MPRedundantFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPRedundantFilter.h index fbab5c18cb517..53c6c611ea1ec 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MPRedundantFilter.h +++ b/L1Trigger/DTTriggerPhase2/interface/MPRedundantFilter.h @@ -27,6 +27,11 @@ class MPRedundantFilter : public MPFilter { const edm::EventSetup& iEventSetup, std::vector& inMPath, std::vector& outMPath) override{}; + void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + std::vector& inSLMPath, + std::vector& inCorMPath, + std::vector& outMPath) override{}; void run(edm::Event& iEvent, const edm::EventSetup& iEventSetup, MuonPathPtrs& inMPath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MPSLFilter.h b/L1Trigger/DTTriggerPhase2/interface/MPSLFilter.h new file mode 100644 index 0000000000000..63d5ec54182c1 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MPSLFilter.h @@ -0,0 +1,71 @@ +#ifndef Phase2L1Trigger_DTTrigger_MPSLFilter_h +#define Phase2L1Trigger_DTTrigger_MPSLFilter_h + +#include "L1Trigger/DTTriggerPhase2/interface/MPFilter.h" +#include "L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h" + +#include +#include + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +// =============================================================================== +// Class declarations +// =============================================================================== + +struct valid_tp_t { + bool valid; + cmsdt::metaPrimitive mp; + valid_tp_t() : valid(false), mp(cmsdt::metaPrimitive()) {} + valid_tp_t(bool valid, cmsdt::metaPrimitive mp) : valid(valid), mp(mp) {} +}; + +using valid_tp_arr_t = std::vector; + +class MPSLFilter : public MPFilter { +public: + // Constructors and destructor + MPSLFilter(const edm::ParameterSet &pset); + ~MPSLFilter() override = default; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPath, + std::vector &outMPath) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inSLMPath, + std::vector &inCorMPath, + std::vector &outMPath) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMPath, + MuonPathPtrs &outMPath) override{}; + + void finish() override; + + // Other public methods + + // Public attributes + void printmP(cmsdt::metaPrimitive mP); + +private: + // Private methods + // std::vector filter(std::map>); + std::vector filter(std::vector mps); + bool isDead(cmsdt::metaPrimitive mp, std::map tps_per_bx); + int killTps(cmsdt::metaPrimitive mp, int bx, std::map &tps_per_bx); + int share_hit(cmsdt::metaPrimitive mp, cmsdt::metaPrimitive mp2); + int match(cmsdt::metaPrimitive mp1, cmsdt::metaPrimitive mp2); + int smaller_chi2(cmsdt::metaPrimitive mp, cmsdt::metaPrimitive mp2); + int get_chi2(cmsdt::metaPrimitive mp); + + // Private attributes + const bool debug_; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPath.h b/L1Trigger/DTTriggerPhase2/interface/MuonPath.h index 9355cf0cbaee6..b23e4d6ef53bb 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MuonPath.h +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPath.h @@ -1,130 +1,134 @@ -#ifndef L1Trigger_DTTriggerPhase2_MuonPath_h -#define L1Trigger_DTTriggerPhase2_MuonPath_h -#include -#include - -#include "L1Trigger/DTTriggerPhase2/interface/DTprimitive.h" - -class MuonPath { -public: - MuonPath(); - MuonPath(DTPrimitivePtrs &ptrPrimitive, int prup = 0, int prdw = 0); - MuonPath(DTPrimitives &ptrPrimitive, int prup = 0, int prdw = 0); - MuonPath(std::shared_ptr &ptr); - virtual ~MuonPath() {} - - // setter methods - void setPrimitive(DTPrimitivePtr &ptr, int layer); - void setNPrimitives(short nprim) { nprimitives_ = nprim; } - void setNPrimitivesUp(short nprim) { nprimitivesUp_ = nprim; } - void setNPrimitivesDown(short nprim) { nprimitivesDown_ = nprim; } - void setCellHorizontalLayout(int layout[4]); - void setCellHorizontalLayout(const int *layout); - void setBaseChannelId(int bch) { baseChannelId_ = bch; } - void setQuality(cmsdt::MP_QUALITY qty) { quality_ = qty; } - void setBxTimeValue(int time); - void setLateralComb(cmsdt::LATERAL_CASES latComb[4]); - void setLateralComb(const cmsdt::LATERAL_CASES *latComb); - void setLateralCombFromPrimitives(void); - - void setHorizPos(float pos) { horizPos_ = pos; } - void setTanPhi(float tanPhi) { tanPhi_ = tanPhi; } - void setChiSquare(float chi) { chiSquare_ = chi; } - void setPhi(float phi) { phi_ = phi; } - void setPhiB(float phib) { phiB_ = phib; } - void setPhiCMSSW(float phi_cmssw) { phicmssw_ = phi_cmssw; } - void setPhiBCMSSW(float phib_cmssw) { phiBcmssw_ = phib_cmssw; } - void setXCoorCell(float x, int cell) { xCoorCell_[cell] = x; } - void setDriftDistance(float dx, int cell) { xDriftDistance_[cell] = dx; } - void setXWirePos(float x, int cell) { xWirePos_[cell] = x; } - void setZWirePos(float z, int cell) { zWirePos_[cell] = z; } - void setTWireTDC(float t, int cell) { tWireTDC_[cell] = t; } - void setRawId(uint32_t id) { rawId_ = id; } - - // getter methods - DTPrimitivePtr primitive(int layer) const { return prim_[layer]; } - short nprimitives() const { return nprimitives_; } - short nprimitivesDown() const { return nprimitivesDown_; } - short nprimitivesUp() const { return nprimitivesUp_; } - const int *cellLayout() const { return cellLayout_; } - int baseChannelId() const { return baseChannelId_; } - cmsdt::MP_QUALITY quality() const { return quality_; } - int bxTimeValue() const { return bxTimeValue_; } - int bxNumId() const { return bxNumId_; } - float tanPhi() const { return tanPhi_; } - const cmsdt::LATERAL_CASES *lateralComb() const { return (lateralComb_); } - float horizPos() const { return horizPos_; } - float chiSquare() const { return chiSquare_; } - float phi() const { return phi_; } - float phiB() const { return phiB_; } - float phi_cmssw() const { return phicmssw_; } - float phiB_cmssw() const { return phiBcmssw_; } - float xCoorCell(int cell) const { return xCoorCell_[cell]; } - float xDriftDistance(int cell) const { return xDriftDistance_[cell]; } - float xWirePos(int cell) const { return xWirePos_[cell]; } - float zWirePos(int cell) const { return zWirePos_[cell]; } - float tWireTDC(int cell) const { return tWireTDC_[cell]; } - uint32_t rawId() const { return rawId_; } - - // Other methods - bool isEqualTo(MuonPath *ptr); - bool isAnalyzable(); - bool completeMP(); - -private: - //------------------------------------------------------------------ - //--- MuonPath's data - //------------------------------------------------------------------ - /* - Primitives that make up the path. The 0th position holds the channel ID of - the lower layer. The order is critical. - */ - DTPrimitivePtrs prim_; //ENSURE that there are no more than 4-8 prims - short nprimitives_; - short nprimitivesUp_; - short nprimitivesDown_; - - /* Horizontal position of each cell (one per layer), in half-cell units, - with respect of the lower layer (layer 0). - */ - int cellLayout_[cmsdt::NUM_LAYERS]; - int baseChannelId_; - - //------------------------------------------------------------------ - //--- Fit results: - //------------------------------------------------------------------ - /* path quality */ - cmsdt::MP_QUALITY quality_; - - /* Lateral combination */ - cmsdt::LATERAL_CASES lateralComb_[cmsdt::NUM_LAYERS]; - - /* BX time value with respect to BX0 of the orbit */ - int bxTimeValue_; - - /* BX number in the orbit */ - int bxNumId_; - - /* Cell parameters */ - float xCoorCell_[cmsdt::NUM_LAYERS_2SL]; // Horizontal position of the hit in each cell - float xDriftDistance_[cmsdt::NUM_LAYERS_2SL]; // Drift distance on the cell (absolute value) - float xWirePos_[cmsdt::NUM_LAYERS_2SL]; - float zWirePos_[cmsdt::NUM_LAYERS_2SL]; - float tWireTDC_[cmsdt::NUM_LAYERS_2SL]; - - float tanPhi_; - float horizPos_; - float chiSquare_; - float phi_; - float phiB_; - float phicmssw_; - float phiBcmssw_; - - uint32_t rawId_; -}; - -typedef std::vector MuonPaths; -typedef std::shared_ptr MuonPathPtr; -typedef std::vector MuonPathPtrs; - -#endif +#ifndef L1Trigger_DTTriggerPhase2_MuonPath_h +#define L1Trigger_DTTriggerPhase2_MuonPath_h +#include +#include + +#include "L1Trigger/DTTriggerPhase2/interface/DTprimitive.h" + +class MuonPath { +public: + MuonPath(); + MuonPath(DTPrimitivePtrs &ptrPrimitive, int prup = 0, int prdw = 0); + MuonPath(DTPrimitives &ptrPrimitive, int prup = 0, int prdw = 0); + MuonPath(std::shared_ptr &ptr); + virtual ~MuonPath() {} + + // setter methods + void setPrimitive(DTPrimitivePtr &ptr, int layer); + void setNPrimitives(short nprim) { nprimitives_ = nprim; } + void setNPrimitivesUp(short nprim) { nprimitivesUp_ = nprim; } + void setNPrimitivesDown(short nprim) { nprimitivesDown_ = nprim; } + void setCellHorizontalLayout(int layout[4]); + void setCellHorizontalLayout(const int *layout); + void setBaseChannelId(int bch) { baseChannelId_ = bch; } + void setMissingLayer(int layer) { missingLayer_ = layer; } + void setQuality(cmsdt::MP_QUALITY qty) { quality_ = qty; } + void setBxTimeValue(int time); + void setLateralComb(cmsdt::LATERAL_CASES latComb[4]); + void setLateralComb(const cmsdt::LATERAL_CASES *latComb); + void setLateralCombFromPrimitives(void); + + void setHorizPos(float pos) { horizPos_ = pos; } + void setTanPhi(float tanPhi) { tanPhi_ = tanPhi; } + void setChiSquare(float chi) { chiSquare_ = chi; } + void setPhi(float phi) { phi_ = phi; } + void setPhiB(float phib) { phiB_ = phib; } + void setPhiCMSSW(float phi_cmssw) { phicmssw_ = phi_cmssw; } + void setPhiBCMSSW(float phib_cmssw) { phiBcmssw_ = phib_cmssw; } + void setXCoorCell(float x, int cell) { xCoorCell_[cell] = x; } + void setDriftDistance(float dx, int cell) { xDriftDistance_[cell] = dx; } + void setXWirePos(float x, int cell) { xWirePos_[cell] = x; } + void setZWirePos(float z, int cell) { zWirePos_[cell] = z; } + void setTWireTDC(float t, int cell) { tWireTDC_[cell] = t; } + void setRawId(uint32_t id) { rawId_ = id; } + + // getter methods + DTPrimitivePtr primitive(int layer) const { return prim_[layer]; } + short nprimitives() const { return nprimitives_; } + short nprimitivesDown() const { return nprimitivesDown_; } + short nprimitivesUp() const { return nprimitivesUp_; } + const int *cellLayout() const { return cellLayout_; } + int baseChannelId() const { return baseChannelId_; } + int missingLayer() const { return missingLayer_; } + cmsdt::MP_QUALITY quality() const { return quality_; } + int bxTimeValue() const { return bxTimeValue_; } + int bxNumId() const { return bxNumId_; } + float tanPhi() const { return tanPhi_; } + const cmsdt::LATERAL_CASES *lateralComb() const { return (lateralComb_); } + float horizPos() const { return horizPos_; } + float chiSquare() const { return chiSquare_; } + float phi() const { return phi_; } + float phiB() const { return phiB_; } + float phi_cmssw() const { return phicmssw_; } + float phiB_cmssw() const { return phiBcmssw_; } + float xCoorCell(int cell) const { return xCoorCell_[cell]; } + float xDriftDistance(int cell) const { return xDriftDistance_[cell]; } + float xWirePos(int cell) const { return xWirePos_[cell]; } + float zWirePos(int cell) const { return zWirePos_[cell]; } + float tWireTDC(int cell) const { return tWireTDC_[cell]; } + uint32_t rawId() const { return rawId_; } + + // Other methods + bool isEqualTo(MuonPath *ptr); + bool isAnalyzable(); + bool completeMP(); + +private: + //------------------------------------------------------------------ + //--- MuonPath's data + //------------------------------------------------------------------ + /* + Primitives that make up the path. The 0th position holds the channel ID of + the lower layer. The order is critical. + */ + DTPrimitivePtrs prim_; //ENSURE that there are no more than 4-8 prims + short nprimitives_; + short nprimitivesUp_; + short nprimitivesDown_; + + /* Horizontal position of each cell (one per layer), in half-cell units, + with respect of the lower layer (layer 0). + */ + int cellLayout_[cmsdt::NUM_LAYERS]; + int baseChannelId_; + + int missingLayer_; + + //------------------------------------------------------------------ + //--- Fit results: + //------------------------------------------------------------------ + /* path quality */ + cmsdt::MP_QUALITY quality_; + + /* Lateral combination */ + cmsdt::LATERAL_CASES lateralComb_[cmsdt::NUM_LAYERS]; + + /* BX time value with respect to BX0 of the orbit */ + int bxTimeValue_; + + /* BX number in the orbit */ + int bxNumId_; + + /* Cell parameters */ + float xCoorCell_[cmsdt::NUM_LAYERS_2SL]; // Horizontal position of the hit in each cell + float xDriftDistance_[cmsdt::NUM_LAYERS_2SL]; // Drift distance on the cell (absolute value) + float xWirePos_[cmsdt::NUM_LAYERS_2SL]; + float zWirePos_[cmsdt::NUM_LAYERS_2SL]; + float tWireTDC_[cmsdt::NUM_LAYERS_2SL]; + + float tanPhi_; + float horizPos_; + float chiSquare_; + float phi_; + float phiB_; + float phicmssw_; + float phiBcmssw_; + + uint32_t rawId_; +}; + +typedef std::vector MuonPaths; +typedef std::shared_ptr MuonPathPtr; +typedef std::vector MuonPathPtrs; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h index 9b4f8dd0e32e3..2096cbfa9ea74 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h @@ -53,6 +53,15 @@ class MuonPathAnalyticAnalyzer : public MuonPathAnalyzer { const edm::EventSetup &iEventSetup, MuonPathPtrs &inMpath, std::vector &metaPrimitives) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities, + std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPaths, + std::vector &outMPaths) override{}; void run(edm::Event &iEvent, const edm::EventSetup &iEventSetup, MuonPathPtrs &inMpath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h index e6cc4e0b58cf0..bc00b779fbcda 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h @@ -19,6 +19,7 @@ #include "L1Trigger/DTTriggerPhase2/interface/MuonPath.h" #include "L1Trigger/DTTriggerPhase2/interface/constants.h" #include "L1Trigger/DTTriggerPhase2/interface/GlobalCoordsObtainer.h" +#include "L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h" #include "Geometry/Records/interface/MuonGeometryRecord.h" #include "Geometry/DTGeometry/interface/DTGeometry.h" @@ -46,6 +47,15 @@ class MuonPathAnalyzer { const edm::EventSetup& iEventSetup, MuonPathPtrs& inMpath, std::vector& metaPrimitives) = 0; + virtual void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + MuonPathPtrs& inMpath, + std::vector& lateralities, + std::vector& metaPrimitives) = 0; + virtual void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + std::vector& inMPaths, + std::vector& outMPaths) = 0; virtual void run(edm::Event& iEvent, const edm::EventSetup& iEventSetup, MuonPathPtrs& inMpath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h index 9f2983b9d0c08..d867437359825 100644 --- a/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h @@ -27,6 +27,15 @@ class MuonPathAnalyzerInChamber : public MuonPathAnalyzer { const edm::EventSetup &iEventSetup, MuonPathPtrs &inMpath, std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities, + std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPaths, + std::vector &outMPaths) override{}; void run(edm::Event &iEvent, const edm::EventSetup &iEventSetup, MuonPathPtrs &inMpath, diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h new file mode 100644 index 0000000000000..c953434be6a71 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h @@ -0,0 +1,72 @@ +#ifndef L1Trigger_DTTriggerPhase2_MuonPathConfirmator_h +#define L1Trigger_DTTriggerPhase2_MuonPathConfirmator_h + +#include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/FrameworkfwdMostUsed.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "L1Trigger/DTTriggerPhase2/interface/constants.h" +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonDetId/interface/DTSuperLayerId.h" +#include "DataFormats/MuonDetId/interface/DTLayerId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/DTDigi/interface/DTDigiCollection.h" + +#include +#include + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +// =============================================================================== +// Class declarations +// =============================================================================== + +class MuonPathConfirmator { +public: + // Constructors and destructor + MuonPathConfirmator(const edm::ParameterSet &pset, edm::ConsumesCollector &iC); + ~MuonPathConfirmator(); + + // Main methods + void initialise(const edm::EventSetup &iEventSetup); + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector inMetaPrimitives, + edm::Handle dtdigis, + std::vector &outMetaPrimitives); + + void finish(); + // Other public methods + +private: + // Private methods + void analyze(cmsdt::metaPrimitive mp, + edm::Handle dtdigis, + std::vector &outMetaPrimitives); + // Private attributes + bool debug_; + double minx_match_2digis_; + //shift + edm::FileInPath shift_filename_; + std::map shiftinfo_; + edm::FileInPath maxdrift_filename_; + int maxdriftinfo_[5][4][14]; + int max_drift_tdc = -1; + + int PARTIALS_PRECISSION = 4; + int SEMICHAMBER_H_PRECISSION = 13 + PARTIALS_PRECISSION; + int SEMICHAMBER_RES_SHR = SEMICHAMBER_H_PRECISSION; + int LYRANDAHALF_RES_SHR = 4; + float SEMICHAMBER_H_REAL = ((235. / 2.) / (16. * 6.5)) * std::pow(2, SEMICHAMBER_H_PRECISSION); + int SEMICHAMBER_H = int(SEMICHAMBER_H_REAL); +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h new file mode 100644 index 0000000000000..fa50ddebae2ca --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h @@ -0,0 +1,74 @@ +#ifndef L1Trigger_DTTriggerPhase2_MuonPathCorFitter_h +#define L1Trigger_DTTriggerPhase2_MuonPathCorFitter_h + +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h" + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +inline bool bxSort(const cmsdt::bx_sl_vector &vA, const cmsdt::bx_sl_vector &vB) { + if (vA.bx > vB.bx) + return true; + else if (vA.bx == vB.bx) + return (vA.sl > vB.sl); + else + return false; +} + +using mp_group = std::vector; + +// =============================================================================== +// Class declarations +// =============================================================================== + +class MuonPathCorFitter : public MuonPathFitter { +public: + // Constructors and destructor + MuonPathCorFitter(const edm::ParameterSet &pset, + edm::ConsumesCollector &iC, + std::shared_ptr &globalcoordsobtainer); + ~MuonPathCorFitter() override; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities, + std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPaths, + std::vector &outMPaths) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + MuonPathPtrs &outMPath) override{}; + + void finish() override; + + // Other public methods + + // luts + edm::FileInPath both_sl_filename_; + +private: + // Private methods + void analyze(mp_group mp, std::vector &metaPrimitives); + void fillLuts(); + int get_rom_addr(mp_group mps, std::vector missing_hits); + bool canCorrelate(cmsdt::metaPrimitive mp_sl1, cmsdt::metaPrimitive mp_sl3); + + // Private attributes + int dT0_correlate_TP_; + + // double chi2Th_; + std::vector> lut_2sl; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h new file mode 100644 index 0000000000000..3c19d831112c7 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h @@ -0,0 +1,115 @@ +#ifndef L1Trigger_DTTriggerPhase2_MuonPathFitter_h +#define L1Trigger_DTTriggerPhase2_MuonPathFitter_h + +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h" +#include "L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h" + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +using coeff_arr_t = std::vector>; +struct coeffs_t { + coeff_arr_t t0; + coeff_arr_t position; + coeff_arr_t slope; + coeffs_t() + : t0(cmsdt::N_COEFFS, std::vector(cmsdt::GENERIC_COEFF_WIDTH, 0)), + position(cmsdt::N_COEFFS, std::vector(cmsdt::GENERIC_COEFF_WIDTH, 0)), + slope(cmsdt::N_COEFFS, std::vector(cmsdt::GENERIC_COEFF_WIDTH, 0)) {} +}; + +struct SLhitP { + int ti; // unsigned(16 downto 0); -- 12 msb = bunch_ctr, 5 lsb = tdc counts, resolution 25/32 ns + int wi; // unsigned(6 downto 0); -- ~ 96 channels per layer + int ly; // unsigned(1 downto 0); -- 4 layers + int wp; // signed(WIREPOS_WIDTH-1 downto 0); +}; + +struct fit_common_in_t { + // int valid; not needed, we will not propagate the mpath to the fitter + std::vector hits; + std::vector hits_valid; // slv(0 to 7) + std::vector lateralities; // slv(0 to 7) + coeffs_t coeffs; + int coarse_bctr; // unsigned(11 downto 0) + int coarse_wirepos; // signed(WIDTH_FULL_POS-1 downto WIREPOS_NORM_LSB_IGNORED); +}; + +struct fit_common_out_t { + int t0; + int slope; + int position; + int chi2; + int valid_fit; + fit_common_out_t() : t0(0), slope(0), position(0), chi2(0), valid_fit(0) {} +}; + +// =============================================================================== +// Class declarations +// =============================================================================== + +class MuonPathFitter : public MuonPathAnalyzer { +public: + // Constructors and destructor + MuonPathFitter(const edm::ParameterSet &pset, + edm::ConsumesCollector &iC, + std::shared_ptr &globalcoordsobtainer); + ~MuonPathFitter() override; + + // Main methods + + // Other public methods + coeffs_t RomDataConvert(std::vector slv, + short COEFF_WIDTH_T0, + short COEFF_WIDTH_POSITION, + short COEFF_WIDTH_SLOPE, + short LOLY, + short HILY); + + bool hasPosRF(int wh, int sec) { return wh > 0 || (wh == 0 && sec % 4 > 1); }; + void setChi2Th(double chi2Th) { chi2Th_ = chi2Th; }; + void setTanPhiTh(double tanPhiTh) { tanPhiTh_ = tanPhiTh; }; + + // Public attributes + DTGeometry const *dtGeo_; + edm::ESGetToken dtGeomH; + + //shift + edm::FileInPath shift_filename_; + std::map shiftinfo_; + + // max drift velocity + edm::FileInPath maxdrift_filename_; + int maxdriftinfo_[5][4][14]; + int max_drift_tdc = -1; + + int get_rom_addr(MuonPathPtr &inMPath, latcomb lats); + fit_common_out_t fit(fit_common_in_t fit_common_in, + int XI_WIDTH, + int COEFF_WIDTH_T0, + int COEFF_WIDTH_POSITION, + int COEFF_WIDTH_SLOPE, + int PRECISSION_T0, + int PRECISSION_POSITION, + int PRECISSION_SLOPE, + int PROD_RESIZE_T0, + int PROD_RESIZE_POSITION, + int PROD_RESIZE_SLOPE, + int MAX_DRIFT_TDC, + int sl); + + double tanPhiTh_; + const bool debug_; + double chi2Th_; + + // global coordinates + std::shared_ptr globalcoordsobtainer_; + +private: + // Private methods + + // Private attributes +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h b/L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h new file mode 100644 index 0000000000000..b15402ff4424e --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h @@ -0,0 +1,69 @@ +#ifndef L1Trigger_DTTriggerPhase2_MuonPathSLFitter_h +#define L1Trigger_DTTriggerPhase2_MuonPathSLFitter_h + +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h" + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +// =============================================================================== +// Class declarations +// =============================================================================== + +class MuonPathSLFitter : public MuonPathFitter { +public: + // Constructors and destructor + MuonPathSLFitter(const edm::ParameterSet &pset, + edm::ConsumesCollector &iC, + std::shared_ptr &globalcoordsobtainer); + ~MuonPathSLFitter() override; + + // Main methods + void initialise(const edm::EventSetup &iEventSetup) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &metaPrimitives) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + std::vector &lateralities, + std::vector &metaPrimitives) override; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPaths, + std::vector &outMPaths) override{}; + void run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &inMpath, + MuonPathPtrs &outMPath) override{}; + + void finish() override; + + // Other public methods + + //shift theta + edm::FileInPath shift_theta_filename_; + std::map shiftthetainfo_; + + // luts + edm::FileInPath sl1_filename_; + edm::FileInPath sl2_filename_; + edm::FileInPath sl3_filename_; + +private: + // Private methods + void analyze(MuonPathPtr &inMPath, lat_vector lat_combs, std::vector &metaPrimitives); + void fillLuts(); + int get_rom_addr(MuonPathPtr &inMPath, latcomb lats); + + // Private attributes + + // double chi2Th_; + std::vector> lut_sl1; + std::vector> lut_sl2; + std::vector> lut_sl3; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h b/L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h new file mode 100644 index 0000000000000..323831af2ee9b --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h @@ -0,0 +1,332 @@ +#ifndef Phase2L1Trigger_DTTrigger_TrapezoidalGrouping_h +#define Phase2L1Trigger_DTTrigger_TrapezoidalGrouping_h + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonDetId/interface/DTSuperLayerId.h" +#include "DataFormats/MuonDetId/interface/DTLayerId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/DTDigi/interface/DTDigiCollection.h" + +#include "L1Trigger/DTTriggerPhase2/interface/MuonPath.h" +#include "L1Trigger/DTTriggerPhase2/interface/constants.h" + +#include "L1Trigger/DTTriggerPhase2/interface/MotherGrouping.h" + +#include +#include +#include + +// =============================================================================== +// Previous definitions and declarations +// =============================================================================== + +/* + Channels are labeled following next schema: + --------------------------------- + | 6 | 7 | 8 | 9 | + --------------------------------- + | 3 | 4 | 5 | + ------------------------- + | 1 | 2 | + ----------------- + | 0 | + --------- +*/ + +inline bool hitWireSort(const DTPrimitive& hit1, const DTPrimitive& hit2) { + int wi1 = hit1.channelId(); + int wi2 = hit2.channelId(); + + if (wi1 < wi2) + return true; + else + return false; +} + +inline bool hitLayerSort(const DTPrimitive& hit1, const DTPrimitive& hit2) { + int lay1 = hit1.layerId(); + int lay2 = hit2.layerId(); + + if (lay1 < lay2) + return true; + else if (lay1 > lay2) + return false; + else + return hitWireSort(hit1, hit2); +} + +inline bool hitTimeSort(const DTPrimitive& hit1, const DTPrimitive& hit2) { + int tdc1 = hit1.tdcTimeStamp(); + int tdc2 = hit2.tdcTimeStamp(); + + if (tdc1 < tdc2) + return true; + else if (tdc1 > tdc2) + return false; + else + return hitLayerSort(hit1, hit2); +} + +namespace dtamgrouping { + /* Cell's combination, following previous labeling, to obtain every possible muon's path. + Others cells combinations imply non straight paths */ + constexpr int CHANNELS_PATH_ARRANGEMENTS[8][4] = { + {0, 1, 3, 6}, {0, 1, 3, 7}, {0, 1, 4, 7}, {0, 1, 4, 8}, {0, 2, 4, 7}, {0, 2, 4, 8}, {0, 2, 5, 8}, {0, 2, 5, 9}}; + + /* For each of the previous cell's combinations, this array stores the associated cell's + displacement, relative to lower layer cell, measured in semi-cell length units */ + + constexpr int CELL_HORIZONTAL_LAYOUTS[8][4] = {{0, -1, -2, -3}, + {0, -1, -2, -1}, + {0, -1, 0, -1}, + {0, -1, 0, 1}, + {0, 1, 0, -1}, + {0, 1, 0, 1}, + {0, 1, 2, 1}, + {0, 1, 2, 3}}; +} // namespace dtamgrouping + +// =============================================================================== +// Class declarations +// =============================================================================== + +class TrapezoidalGrouping : public MotherGrouping { +public: + // Constructors and destructor + TrapezoidalGrouping(const edm::ParameterSet& pset, edm::ConsumesCollector& iC); + ~TrapezoidalGrouping() override; + + // Main methods + void initialise(const edm::EventSetup& iEventSetup) override; + void run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + const DTDigiCollection& digis, + MuonPathPtrs& outMpath) override; + void finish() override; + + // Other public methods + + // Public attributes + +private: + // Private methods + void setInChannels(const DTDigiCollection* digi, int sl); + std::vector group_hits(DTPrimitive pivot_hit, + std::vector input_paths, + DTPrimitives hits_per_cell, + DTPrimitives& hits_in_trapezoid); + + // Private attributes + const bool debug_; + + DTPrimitives muxInChannels_[cmsdt::NUM_CELLS_PER_BLOCK]; + DTPrimitives channelIn_[cmsdt::NUM_LAYERS][cmsdt::NUM_CH_PER_LAYER]; + DTPrimitives all_hits; + DTPrimitives chInDummy_; + int prevTDCTimeStamps_[4]; + int currentBaseChannel_; + + // The trapezoid is as follows: + // [ 0 ][ 1 ][ 2 ][ 3 ][ 4 ][ 5 ][ 6 ][ 7 ][ 8 ] + + // And maps to the physical cells as follows: + + // Pivot in layer 1 = "00" + // [ 5 ][ 6 ][ 7 ][ 8 ] Layer C + // [ 2 ][ 3 ][ 4 ] Layer B + // [ 0 ][ 1 ] Layer A + // Pivot + + // Pivot in layer 2 = "01" + // [ 2 ][ 3 ][ 4 ] Layer B + // [ 0 ][ 1 ] Layer A + // Pivot + // [ 6,8 ][ 5,7 ] Layer C + + // Pivot in layer 3 = "10" + // [ 6,8 ][ 5,7 ] Layer C + // Pivot + // [ 0 ][ 1 ] Layer A + // [ 2 ][ 3 ][ 4 ] Layer B + + // Pivot in layer 4 = "11" + // Pivot + // [ 0 ][ 1 ] Layer A + // [ 2 ][ 3 ][ 4 ] Layer B + // [ 5 ][ 6 ][ 7 ][ 8 ] Layer C + + short trapezoid_vertical_mapping[4][9] = {{1, 1, 2, 2, 2, 3, 3, 3, 3}, + {1, 1, 2, 2, 2, -1, -1, -1, -1}, + {-1, -1, -2, -2, -2, 1, 1, 1, 1}, + {-1, -1, -2, -2, -2, -3, -3, -3, -3}}; + + short trapezoid_horizontal_mapping[4][9] = {{0, 1, -1, 0, 1, -1, 0, 1, 2}, + {-1, 0, -1, 0, 1, 0, -1, 0, -1}, + {0, 1, -1, 0, 1, 1, 0, 1, 0}, + {-1, 0, -1, 0, 1, -2, -1, 0, 1}}; + + // Task list + // 4 hit candidates + // 0 => (0,2,5), 1 => (0,2,6), 2 => (0,3,6), 3 => (0,3,7), + // 4 => (1,3,6), 5 => (1,3,7), 6 => (1,4,7), 7 => (1,4,8), + // the rest are 3-hit candidates, last value not used + // 8 => (0,2,0), 9 => (0,3,0), 10 => (1,3,0), 11 => (1,4,0), + // 12 => (0,5,0), 13 => (0,6,0), 14 => (0,7,0), 15 => (1,6,0), + // 16 => (1,7,0), 17 => (1,8,0), 18 => (2,5,0), 19 => (2,6,0), + // 20 => (3,6,0), 21 => (3,7,0), 22 => (4,7,0), 23 => (4,8,0) + + std::vector> task_list = {// 4-hit + {0, 2, 5}, + {0, 2, 6}, + {0, 3, 6}, + {0, 3, 7}, + {1, 3, 6}, + {1, 3, 7}, + {1, 4, 7}, + {1, 4, 8}, + // 3-hit + {0, 2}, + {0, 3}, + {1, 3}, + {1, 4}, + {0, 5}, + {0, 6}, + {0, 7}, + {1, 6}, + {1, 7}, + {1, 8}, + {2, 5}, + {2, 6}, + {3, 6}, + {3, 7}, + {4, 7}, + {4, 8}}; + + int CELL_HORIZONTAL_LAYOUTS_PER_TASK[4][24][4] = { // pivoting over layer 1 + {// all layers available + {0, 0, 0, -1}, + {0, 0, 1, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {1, 0, 0, -1}, + {1, 0, 1, -1}, + {1, 1, 0, -1}, + {1, 1, 1, -1}, + // layer 4 missing + {0, 0, 0, -1}, + {0, 1, 0, -1}, + {1, 0, 0, -1}, + {1, 1, 0, -1}, + // layer 3 missing + {0, 0, 0, -1}, + {0, 0, 1, -1}, + {0, 1, 1, -1}, + {1, 0, 0, -1}, + {1, 0, 1, -1}, + {1, 1, 1, -1}, + // layer 2 missing + {0, 0, 0, -1}, + {0, 0, 1, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {1, 1, 0, -1}, + {1, 1, 1, -1}}, + // pivoting over layer 2 + {// all layers available + {0, 0, 0, -1}, + {1, 0, 0, -1}, + {1, 0, 1, -1}, + {0, 0, 1, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {1, 1, 1, -1}, + // layer 1 missing + {0, 0, 0, -1}, + {0, 0, 1, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + // layer 4 missing + {0, 0, 0, -1}, + {1, 0, 0, -1}, + {0, 0, 0, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {1, 1, 0, -1}, + // layer 3 missing + {0, 0, 0, -1}, + {1, 0, 0, -1}, + {1, 0, 1, -1}, + {0, 0, 1, -1}, + {0, 1, 1, -1}, + {1, 1, 1, -1}}, + // pivoting over layer 3 + {// all layers available + {1, 1, 1, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {1, 0, 0, -1}, + {1, 0, 1, -1}, + {0, 0, 1, -1}, + {0, 0, 0, -1}, + // layer 4 missing + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {1, 0, 0, -1}, + {0, 0, 0, -1}, + // layer 1 missing + {0, 1, 1, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {0, 0, 0, -1}, + {0, 0, 1, -1}, + {0, 0, 0, -1}, + // layer 2 missing + {1, 1, 1, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {0, 1, 1, -1}, + {0, 0, 1, -1}, + {0, 0, 0, -1}}, + // pivoting over layer 4 + {// all layers available + {1, 1, 1, -1}, + {0, 1, 1, -1}, + {1, 0, 1, -1}, + {0, 0, 1, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {1, 0, 0, -1}, + {0, 0, 0, -1}, + // layer 1 missing + {0, 1, 1, -1}, + {0, 0, 1, -1}, + {0, 1, 0, -1}, + {0, 0, 0, -1}, + // layer 2 missing + {1, 1, 1, -1}, + {0, 1, 1, -1}, + {0, 0, 1, -1}, + {1, 1, 0, -1}, + {0, 1, 0, -1}, + {0, 0, 0, -1}, + // layer 3 missing + {1, 1, 1, -1}, + {0, 1, 1, -1}, + {1, 0, 1, -1}, + {0, 0, 1, -1}, + {1, 0, 0, -1}, + {0, 0, 0, -1}}}; + + int MISSING_LAYER_LAYOUTS_PER_TASK[4][24] = { + {-1, -1, -1, -1, -1, -1, -1, -1, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1}, + {-1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2}, + {-1, -1, -1, -1, -1, -1, -1, -1, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1}, + {-1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2}}; +}; + +#endif diff --git a/L1Trigger/DTTriggerPhase2/interface/constants.h b/L1Trigger/DTTriggerPhase2/interface/constants.h index ac69c152d74fe..454fd0abcac63 100644 --- a/L1Trigger/DTTriggerPhase2/interface/constants.h +++ b/L1Trigger/DTTriggerPhase2/interface/constants.h @@ -20,6 +20,8 @@ #ifndef L1Trigger_DTTriggerPhase2_constants_h #define L1Trigger_DTTriggerPhase2_constants_h #include +#include +#include // Compiler option to select program mode: PRUEBA_MEZCLADOR, PRUEBA_ANALIZADOR, // or NONE @@ -119,7 +121,43 @@ namespace cmsdt { lat8(l8), index(idx), rpcFlag(rpc) {} - + metaPrimitive() + : rawId(0), + t0(0), + x(0), + tanPhi(0), + phi(0), + phiB(0), + phi_cmssw(0), + phiB_cmssw(0), + chi2(0), + quality(0), + wi1(0), + tdc1(0), + lat1(0), + wi2(0), + tdc2(0), + lat2(0), + wi3(0), + tdc3(0), + lat3(0), + wi4(0), + tdc4(0), + lat4(0), + wi5(0), + tdc5(0), + lat5(0), + wi6(0), + tdc6(0), + lat6(0), + wi7(0), + tdc7(0), + lat7(0), + wi8(0), + tdc8(0), + lat8(0), + index(0), + rpcFlag(0) {} uint32_t rawId; double t0; double x; @@ -157,10 +195,12 @@ namespace cmsdt { int index; int rpcFlag = 0; }; + struct PARTIAL_LATQ_TYPE { bool latQValid; int bxValue; }; + struct LATQ_TYPE { bool valid; int bxValue; @@ -168,6 +208,12 @@ namespace cmsdt { MP_QUALITY quality; }; + struct bx_sl_vector { + int bx; + std::vector mps; + int sl; + }; + enum algo { Standard = 0, PseudoBayes = 1, HoughTrans = 2 }; enum scenario { MC = 0, DATA = 1, SLICE_TEST = 2 }; @@ -175,11 +221,92 @@ namespace cmsdt { /* En nanosegundos */ constexpr int LHC_CLK_FREQ = 25; + /* mixer constants */ + // Hits can be separated up to 9 frames, with 2 BXs per frame + // | | | | | | | | | | | | | | | | | | + // F1 F2 F3 F4 F5 F6 F7 F8 F9 + constexpr int BX_PER_FRAME = 2; + constexpr int MAX_FRAME_DIF = 8; + constexpr int PATHFINDER_INPUT_HITS_LIMIT = 8; + + /* laterality provider */ + constexpr int LAT_TOTAL_BITS = 9; // tdc counts from 0 to 512 + constexpr int LAT_MSB_BITS = 6; + constexpr int TIME_TO_TDC_COUNTS = 32; + + constexpr int LAT_P0_4H = 1; + constexpr int LAT_P1_4H = 31; + constexpr int LAT_P2_4H = 40; + + constexpr int LAT_P0_3H = 24; + constexpr int LAT_P1_3H = 27; + constexpr int LAT_P2_3H = 30; + + /* Fitting */ + + constexpr int SL1_CELLS_OFFSET = 48; + + constexpr int N_COEFFS = 8; + constexpr int GENERIC_COEFF_WIDTH = 20; + constexpr int WIDTH_FULL_TIME = 17; + constexpr int WIDTH_COARSED_TIME = 12; + constexpr int WIDTH_DIFBX = 5; + constexpr int WIDTH_FULL_POS = 17; + constexpr int WIDTH_FULL_SLOPE = 14; + constexpr int WIDTH_FULL_CHI2 = 16; + constexpr int WIREPOS_WIDTH = 17; + constexpr int WIREPOS_NORM_LSB_IGNORED = 9; + constexpr int WIDTH_POS_SLOPE_CORR = 9; + + constexpr int XI_SL_WIDTH = 12; + + constexpr int COEFF_WIDTH_SL_T0 = 15; + constexpr int COEFF_WIDTH_SL_POSITION = 18; + constexpr int COEFF_WIDTH_SL2_POSITION = 15; + constexpr int COEFF_WIDTH_SL_SLOPE = 18; + + constexpr int PRECISSION_SL_T0 = 13; + constexpr int PRECISSION_SL_POSITION = 13; + constexpr int PRECISSION_SL_SLOPE = 13; + + constexpr int PROD_RESIZE_SL_T0 = 28; + constexpr int PROD_RESIZE_SL_POSITION = 30; + constexpr int PROD_RESIZE_SL2_POSITION = 27; + constexpr int PROD_RESIZE_SL_SLOPE = 30; + + constexpr int XI_COR_WIDTH = 14; + + constexpr int COEFF_WIDTH_COR_T0 = 15; + constexpr int COEFF_WIDTH_COR_POSITION = 15; + constexpr int COEFF_WIDTH_COR_SLOPE = 15; + + constexpr int PRECISSION_COR_T0 = 15; + constexpr int PRECISSION_COR_POSITION = 15; + constexpr int PRECISSION_COR_SLOPE = 15; + + constexpr int PROD_RESIZE_COR_T0 = 30; + constexpr int PROD_RESIZE_COR_POSITION = 30; + constexpr int PROD_RESIZE_COR_SLOPE = 29; + + constexpr int T0_CUT_TOLERANCE = 0; + + // Filtering + constexpr int FSEG_T0_BX_LSB = 2; + constexpr int FSEG_T0_DISCARD_LSB = 5; + constexpr int FSEG_T0_SIZE = FSEG_T0_BX_LSB + (5 - FSEG_T0_DISCARD_LSB); + constexpr int FSEG_POS_DISCARD_LSB = 9; + constexpr int FSEG_POS_SIZE = WIDTH_FULL_POS - FSEG_POS_DISCARD_LSB; + constexpr int FSEG_SLOPE_DISCARD_LSB = 9; + constexpr int FSEG_SLOPE_SIZE = WIDTH_FULL_SLOPE - FSEG_SLOPE_DISCARD_LSB; + constexpr int SLFILT_MAX_SEG1T0_TO_SEG2ARRIVAL = 24; + /* Adimensional */ constexpr int MAX_BX_IDX = 3564; // In ns (maximum drift time inside the cell) constexpr float MAXDRIFT = 387; + constexpr float MAXDRIFTTDC = 496; // we could make this value depend on the chamber, to be seen + // In mm (cell dimmensions) constexpr int CELL_HEIGHT = 13; constexpr float CELL_SEMIHEIGHT = 6.5; @@ -190,6 +317,9 @@ namespace cmsdt { // With 4 bits for the decimal part constexpr int DRIFT_SPEED_X4 = 889; // 55.5 * 2 ** 4 + // slope conversion 1 LSB = (v_drift) x (1 tdc count) / (1 semicell_h * 16) ~= 0.4e-3 + constexpr float SLOPE_LSB = ((float)CELL_SEMILENGTH / MAXDRIFTTDC) * (1) / (CELL_SEMIHEIGHT * 16.); + // distance between SLs, cm constexpr float VERT_PHI1_PHI3 = 23.5; @@ -199,6 +329,15 @@ namespace cmsdt { // distance between center of the chamber and each SL in mm, 2 bit precision for the decimal part constexpr int CH_CENTER_TO_MID_SL_X2 = 470; // 117.5 * 2 ** 2 + // max difference in BX to even try to correlate + constexpr int MAX_BX_FOR_COR = 2; + + // max number of TPs to store per BX + constexpr int MAX_PRIM_PER_BX_FOR_COR = 6; + + // max number of TPs to correlate and perform the refitting + constexpr int MAX_PRIM_FOR_COR = 12; + /* This is the maximum value than internal time can take. This is because internal time is cyclical due to the limited size of the time counters and diff --git a/L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h b/L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h new file mode 100644 index 0000000000000..f612f773e00ff --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h @@ -0,0 +1,19 @@ +#ifndef L1Trigger_DTTriggerPhase2_vhdl_h +#define L1Trigger_DTTriggerPhase2_vhdl_h + +#include +#include +#include + +// "à la vhdl" functions +std::vector vhdl_slice(std::vector v, int upper, int lower); +int vhdl_unsigned_to_int(std::vector v); +int vhdl_signed_to_int(std::vector v); +void vhdl_int_to_unsigned(int value, std::vector &v); +void vhdl_int_to_signed(int value, std::vector &v); +void vhdl_resize_unsigned(std::vector &v, int new_size); +void vhdl_resize_signed(std::vector &v, int new_size); +bool vhdl_resize_signed_ok(std::vector v, int new_size); +bool vhdl_resize_unsigned_ok(std::vector v, int new_size); + +#endif diff --git a/L1Trigger/DTTriggerPhase2/plugins/DTTrigPhase2Prod.cc b/L1Trigger/DTTriggerPhase2/plugins/DTTrigPhase2Prod.cc index bd7029041f1ae..608528ed88818 100644 --- a/L1Trigger/DTTriggerPhase2/plugins/DTTrigPhase2Prod.cc +++ b/L1Trigger/DTTriggerPhase2/plugins/DTTrigPhase2Prod.cc @@ -23,14 +23,22 @@ #include "L1Trigger/DTTriggerPhase2/interface/constants.h" #include "L1Trigger/DTTriggerPhase2/interface/MotherGrouping.h" -#include "L1Trigger/DTTriggerPhase2/interface/InitialGrouping.h" +#include "L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h" #include "L1Trigger/DTTriggerPhase2/interface/HoughGrouping.h" #include "L1Trigger/DTTriggerPhase2/interface/PseudoBayesGrouping.h" +#include "L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h" +#include "L1Trigger/DTTriggerPhase2/interface/LateralityBasicProvider.h" +#include "L1Trigger/DTTriggerPhase2/interface/LateralityCoarsedProvider.h" #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzer.h" +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h" +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h" #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyticAnalyzer.h" #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAnalyzerInChamber.h" #include "L1Trigger/DTTriggerPhase2/interface/MuonPathAssociator.h" +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h" #include "L1Trigger/DTTriggerPhase2/interface/MPFilter.h" +#include "L1Trigger/DTTriggerPhase2/interface/MPSLFilter.h" +#include "L1Trigger/DTTriggerPhase2/interface/MPCorFilter.h" #include "L1Trigger/DTTriggerPhase2/interface/MPQualityEnhancerFilter.h" #include "L1Trigger/DTTriggerPhase2/interface/MPRedundantFilter.h" #include "L1Trigger/DTTriggerPhase2/interface/MPCleanHitsFilter.h" @@ -95,7 +103,9 @@ class DTTrigPhase2Prod : public edm::stream::EDProducer<> { bool outer(const metaPrimitive& mp) const; bool inner(const metaPrimitive& mp) const; void printmP(const std::string& ss, const metaPrimitive& mP) const; + void printmP(const metaPrimitive& mP) const; void printmPC(const std::string& ss, const metaPrimitive& mP) const; + void printmPC(const metaPrimitive& mP) const; bool hasPosRF(int wh, int sec) const; // Getter-methods @@ -127,6 +137,15 @@ class DTTrigPhase2Prod : public edm::stream::EDProducer<> { int df_extended_; int max_index_; + bool output_mixer_; + bool output_latpredictor_; + bool output_slfitter_; + bool output_slfilter_; + bool output_confirmed_; + bool output_matcher_; + bool skip_processing_; + bool allow_confirmation_; + // ParameterSet edm::EDGetTokenT dtDigisToken_; edm::EDGetTokenT rpcRecHitsLabel_; @@ -135,11 +154,14 @@ class DTTrigPhase2Prod : public edm::stream::EDProducer<> { int algo_; // Grouping code std::unique_ptr grouping_obj_; std::unique_ptr mpathanalyzer_; + std::unique_ptr latprovider_; std::unique_ptr mpathqualityenhancer_; std::unique_ptr mpathqualityenhancerbayes_; std::unique_ptr mpathredundantfilter_; std::unique_ptr mpathhitsfilter_; - std::unique_ptr mpathassociator_; + std::unique_ptr mpathassociator_; + std::unique_ptr mpathconfirmator_; + std::unique_ptr mpathcorfilter_; std::shared_ptr globalcoordsobtainer_; // Buffering @@ -192,6 +214,16 @@ DTTrigPhase2Prod::DTTrigPhase2Prod(const ParameterSet& pset) // Choosing grouping scheme: algo_ = pset.getParameter("algo"); + // shortcuts + + output_mixer_ = pset.getParameter("output_mixer"); + output_latpredictor_ = pset.getParameter("output_latpredictor"); + output_slfitter_ = pset.getParameter("output_slfitter"); + output_slfilter_ = pset.getParameter("output_slfilter"); + output_confirmed_ = pset.getParameter("output_confirmed"); + output_matcher_ = pset.getParameter("output_matcher"); + allow_confirmation_ = pset.getParameter("allow_confirmation"); + edm::ConsumesCollector consumesColl(consumesCollector()); globalcoordsobtainer_ = std::make_shared(pset); globalcoordsobtainer_->generate_luts(); @@ -203,13 +235,14 @@ DTTrigPhase2Prod::DTTrigPhase2Prod(const ParameterSet& pset) grouping_obj_ = std::make_unique(pset.getParameter("HoughGrouping"), consumesColl); } else { - grouping_obj_ = std::make_unique(pset, consumesColl); + grouping_obj_ = std::make_unique(pset, consumesColl); } if (algo_ == Standard) { if (debug_) LogDebug("DTTrigPhase2Prod") << "DTp2:constructor: JM analyzer"; - mpathanalyzer_ = std::make_unique(pset, consumesColl, globalcoordsobtainer_); + mpathanalyzer_ = std::make_unique(pset, consumesColl, globalcoordsobtainer_); + latprovider_ = std::make_unique(pset, consumesColl); } else { if (debug_) LogDebug("DTTrigPhase2Prod") << "DTp2:constructor: Full chamber analyzer"; @@ -221,11 +254,13 @@ DTTrigPhase2Prod::DTTrigPhase2Prod(const ParameterSet& pset) superCellhalfspacewidth_ = pset.getParameter("superCellspacewidth") / 2; superCelltimewidth_ = pset.getParameter("superCelltimewidth"); - mpathqualityenhancer_ = std::make_unique(pset); + mpathqualityenhancer_ = std::make_unique(pset); mpathqualityenhancerbayes_ = std::make_unique(pset); mpathredundantfilter_ = std::make_unique(pset); mpathhitsfilter_ = std::make_unique(pset); - mpathassociator_ = std::make_unique(pset, consumesColl, globalcoordsobtainer_); + mpathconfirmator_ = std::make_unique(pset, consumesColl); + mpathassociator_ = std::make_unique(pset, consumesColl, globalcoordsobtainer_); + mpathcorfilter_ = std::make_unique(pset); rpc_integrator_ = std::make_unique(pset, consumesColl); dtGeomH = esConsumes(); @@ -249,6 +284,7 @@ void DTTrigPhase2Prod::beginRun(edm::Run const& iRun, const edm::EventSetup& iEv mpathqualityenhancerbayes_->initialise(iEventSetup); // Filter object initialisation mpathhitsfilter_->initialise(iEventSetup); mpathassociator_->initialise(iEventSetup); // Associator object initialisation + mpathcorfilter_->initialise(iEventSetup); if (auto geom = iEventSetup.getHandle(dtGeomH)) { dtGeo_ = &(*geom); @@ -285,7 +321,7 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { else if (debug_) LogDebug("DTTrigPhase2Prod") << "produce - Getting and grouping digis per chamber."; - MuonPathPtrs muonpaths; + std::map muonpaths; for (const auto& ich : dtGeo_->chambers()) { // The code inside this for loop would ideally later fit inside a trigger unit (in principle, a DT station) of the future Phase 2 DT Trigger. const DTChamber* chamb = ich; @@ -328,88 +364,314 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { // the groupings access it, it's not needed to "collect" the final products). while (!superCells.empty()) { - grouping_obj_->run(iEvent, iEventSetup, *(superCells.back()), muonpaths); + grouping_obj_->run(iEvent, iEventSetup, *(superCells.back()), muonpaths[chid.rawId()]); superCells.pop_back(); } } else { - grouping_obj_->run(iEvent, iEventSetup, (*dmit).second, muonpaths); + grouping_obj_->run(iEvent, iEventSetup, (*dmit).second, muonpaths[chid.rawId()]); } } digiMap.clear(); if (dump_) { - for (unsigned int i = 0; i < muonpaths.size(); i++) { - stringstream ss; - ss << iEvent.id().event() << " mpath " << i << ": "; - for (int lay = 0; lay < muonpaths.at(i)->nprimitives(); lay++) - ss << muonpaths.at(i)->primitive(lay)->channelId() << " "; - for (int lay = 0; lay < muonpaths.at(i)->nprimitives(); lay++) - ss << muonpaths.at(i)->primitive(lay)->tdcTimeStamp() << " "; - for (int lay = 0; lay < muonpaths.at(i)->nprimitives(); lay++) - ss << muonpaths.at(i)->primitive(lay)->laterality() << " "; - LogInfo("DTTrigPhase2Prod") << ss.str(); + for (auto& ch_muonpaths : muonpaths) { + for (unsigned int i = 0; i < ch_muonpaths.second.size(); i++) { + stringstream ss; + ss << iEvent.id().event() << " mpath " << i << ": "; + for (int lay = 0; lay < ch_muonpaths.second.at(i)->nprimitives(); lay++) + ss << ch_muonpaths.second.at(i)->primitive(lay)->channelId() << " "; + for (int lay = 0; lay < ch_muonpaths.second.at(i)->nprimitives(); lay++) + ss << ch_muonpaths.second.at(i)->primitive(lay)->tdcTimeStamp() << " "; + for (int lay = 0; lay < ch_muonpaths.second.at(i)->nprimitives(); lay++) + ss << ch_muonpaths.second.at(i)->primitive(lay)->laterality() << " "; + LogInfo("DTTrigPhase2Prod") << ss.str(); + } + } + } + + std::map> lateralities; + if (!output_mixer_) { + for (auto& ch_muonpaths : muonpaths) { + if (algo_ == Standard) { + latprovider_->run(iEvent, iEventSetup, ch_muonpaths.second, lateralities[ch_muonpaths.first]); + } } } // FILTER GROUPING - MuonPathPtrs filteredmuonpaths; - if (algo_ == Standard) { - mpathredundantfilter_->run(iEvent, iEventSetup, muonpaths, filteredmuonpaths); - } else { - mpathhitsfilter_->run(iEvent, iEventSetup, muonpaths, filteredmuonpaths); + std::map filteredmuonpaths; + for (auto& ch_muonpaths : muonpaths) { + if (algo_ == Standard) { + mpathredundantfilter_->run(iEvent, iEventSetup, ch_muonpaths.second, filteredmuonpaths[ch_muonpaths.first]); + } else { + mpathhitsfilter_->run(iEvent, iEventSetup, ch_muonpaths.second, filteredmuonpaths[ch_muonpaths.first]); + } } if (dump_) { - for (unsigned int i = 0; i < filteredmuonpaths.size(); i++) { - stringstream ss; - ss << iEvent.id().event() << " filt. mpath " << i << ": "; - for (int lay = 0; lay < filteredmuonpaths.at(i)->nprimitives(); lay++) - ss << filteredmuonpaths.at(i)->primitive(lay)->channelId() << " "; - for (int lay = 0; lay < filteredmuonpaths.at(i)->nprimitives(); lay++) - ss << filteredmuonpaths.at(i)->primitive(lay)->tdcTimeStamp() << " "; - LogInfo("DTTrigPhase2Prod") << ss.str(); + for (auto& ch_filteredmuonpaths : filteredmuonpaths) { + for (unsigned int i = 0; i < ch_filteredmuonpaths.second.size(); i++) { + stringstream ss; + ss << iEvent.id().event() << " filt. mpath " << i << ": "; + for (int lay = 0; lay < ch_filteredmuonpaths.second.at(i)->nprimitives(); lay++) + ss << ch_filteredmuonpaths.second.at(i)->primitive(lay)->channelId() << " "; + for (int lay = 0; lay < ch_filteredmuonpaths.second.at(i)->nprimitives(); lay++) + ss << ch_filteredmuonpaths.second.at(i)->primitive(lay)->tdcTimeStamp() << " "; + LogInfo("DTTrigPhase2Prod") << ss.str(); + } } } + skip_processing_ = output_mixer_ || output_latpredictor_; + /////////////////////////////////////////// /// Fitting SECTION; /////////////////////////////////////////// - if (debug_) - LogDebug("DTTrigPhase2Prod") << "MUON PATHS found: " << muonpaths.size() << " (" << filteredmuonpaths.size() - << ") in event " << iEvent.id().event(); + if (debug_) { + for (auto& ch_muonpaths : muonpaths) { + LogDebug("DTTrigPhase2Prod") << "MUON PATHS found: " << ch_muonpaths.second.size() << " (" + << filteredmuonpaths[ch_muonpaths.first].size() << ") in event " + << iEvent.id().event(); + } + } if (debug_) LogDebug("DTTrigPhase2Prod") << "filling NmetaPrimtives" << std::endl; - std::vector metaPrimitives; - MuonPathPtrs outmpaths; + std::map> metaPrimitives; + std::map outmpaths; if (algo_ == Standard) { if (debug_) LogDebug("DTTrigPhase2Prod") << "Fitting 1SL "; - mpathanalyzer_->run(iEvent, iEventSetup, filteredmuonpaths, metaPrimitives); + for (auto& ch_muonpaths : muonpaths) { // FIXME, do we need filtered muonpaths? + if (!output_mixer_ && !output_latpredictor_) + mpathanalyzer_->run(iEvent, + iEventSetup, + ch_muonpaths.second, + lateralities[ch_muonpaths.first], + metaPrimitives[ch_muonpaths.first]); + else if (output_mixer_) { + for (auto& inMPath : ch_muonpaths.second) { + auto sl = inMPath->primitive(0)->superLayerId(); // 0, 1, 2 + int selected_lay = 1; + if (inMPath->primitive(0)->tdcTimeStamp() != -1) + selected_lay = 0; + int dumLayId = inMPath->primitive(selected_lay)->cameraId(); + auto dtDumlayerId = DTLayerId(dumLayId); + DTSuperLayerId MuonPathSLId(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), sl + 1); + if (sl == 0) + metaPrimitives[ch_muonpaths.first].emplace_back(metaPrimitive({MuonPathSLId.rawId(), + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + -1, + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + -1, + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + -1, + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1})); + else + metaPrimitives[ch_muonpaths.first].emplace_back(metaPrimitive({MuonPathSLId.rawId(), + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + -1, + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + -1, + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + -1, + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + -1, + -1})); + } + } else if (output_latpredictor_) { + int imp = -1; + for (auto& inMPath : ch_muonpaths.second) { + imp++; + auto sl = inMPath->primitive(0)->superLayerId(); // 0, 1, 2 + int selected_lay = 1; + if (inMPath->primitive(0)->tdcTimeStamp() != -1) + selected_lay = 0; + int dumLayId = inMPath->primitive(selected_lay)->cameraId(); + auto dtDumlayerId = DTLayerId(dumLayId); + DTSuperLayerId MuonPathSLId(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), sl + 1); + for (auto& latcomb : lateralities[ch_muonpaths.first][imp]) { + if (sl == 0) + metaPrimitives[ch_muonpaths.first].emplace_back(metaPrimitive({MuonPathSLId.rawId(), + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + latcomb[0], + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + latcomb[1], + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + latcomb[2], + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + latcomb[3], + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1})); + else + metaPrimitives[ch_muonpaths.first].emplace_back(metaPrimitive({MuonPathSLId.rawId(), + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + latcomb[0], + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + latcomb[1], + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + latcomb[2], + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + latcomb[3], + -1})); + } + } + } + } } else { // implementation for advanced (2SL) grouping, no filter required.. if (debug_) LogDebug("DTTrigPhase2Prod") << "Fitting 2SL at once "; - mpathanalyzer_->run(iEvent, iEventSetup, muonpaths, outmpaths); + for (auto& ch_muonpaths : muonpaths) { + mpathanalyzer_->run(iEvent, iEventSetup, ch_muonpaths.second, outmpaths[ch_muonpaths.first]); + } } + skip_processing_ = skip_processing_ || output_slfitter_; + if (dump_) { - for (unsigned int i = 0; i < outmpaths.size(); i++) { - LogInfo("DTTrigPhase2Prod") << iEvent.id().event() << " mp " << i << ": " << outmpaths.at(i)->bxTimeValue() << " " - << outmpaths.at(i)->horizPos() << " " << outmpaths.at(i)->tanPhi() << " " - << outmpaths.at(i)->phi() << " " << outmpaths.at(i)->phiB() << " " - << outmpaths.at(i)->quality() << " " << outmpaths.at(i)->chiSquare(); + for (auto& ch_outmpaths : outmpaths) { + for (unsigned int i = 0; i < ch_outmpaths.second.size(); i++) { + LogInfo("DTTrigPhase2Prod") << iEvent.id().event() << " mp " << i << ": " + << ch_outmpaths.second.at(i)->bxTimeValue() << " " + << ch_outmpaths.second.at(i)->horizPos() << " " + << ch_outmpaths.second.at(i)->tanPhi() << " " << ch_outmpaths.second.at(i)->phi() + << " " << ch_outmpaths.second.at(i)->phiB() << " " + << ch_outmpaths.second.at(i)->quality() << " " + << ch_outmpaths.second.at(i)->chiSquare(); + } } - for (unsigned int i = 0; i < metaPrimitives.size(); i++) { - stringstream ss; - ss << iEvent.id().event() << " mp " << i << ": "; - printmP(ss.str(), metaPrimitives.at(i)); + for (auto& ch_metaPrimitives : metaPrimitives) { + for (unsigned int i = 0; i < ch_metaPrimitives.second.size(); i++) { + stringstream ss; + ss << iEvent.id().event() << " mp " << i << ": "; + printmP(ss.str(), ch_metaPrimitives.second.at(i)); + } } } muonpaths.clear(); filteredmuonpaths.clear(); + ///////////////////////////////////// + //// CONFIRMATION: + ///////////////////////////////////// + + std::map> confirmedMetaPrimitives; + for (auto& ch_metaPrimitives : metaPrimitives) { + if (!skip_processing_ && allow_confirmation_) + mpathconfirmator_->run( + iEvent, iEventSetup, ch_metaPrimitives.second, dtdigis, confirmedMetaPrimitives[ch_metaPrimitives.first]); + else + for (auto& mp : ch_metaPrimitives.second) { + confirmedMetaPrimitives[ch_metaPrimitives.first].push_back(mp); + } + } + + metaPrimitives.clear(); + skip_processing_ = skip_processing_ || output_confirmed_; + ///////////////////////////////////// // FILTER SECTIONS: //////////////////////////////////// @@ -417,25 +679,38 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { if (debug_) LogDebug("DTTrigPhase2Prod") << "declaring new vector for filtered" << std::endl; - std::vector filteredMetaPrimitives; + std::map> filteredMetaPrimitives; if (algo_ == Standard) - mpathqualityenhancer_->run(iEvent, iEventSetup, metaPrimitives, filteredMetaPrimitives); - + for (auto& ch_confirmedMetaPrimitives : confirmedMetaPrimitives) { + if (!skip_processing_) + mpathqualityenhancer_->run(iEvent, + iEventSetup, + ch_confirmedMetaPrimitives.second, + filteredMetaPrimitives[ch_confirmedMetaPrimitives.first]); + else + for (auto& mp : ch_confirmedMetaPrimitives.second) { + filteredMetaPrimitives[ch_confirmedMetaPrimitives.first].push_back(mp); + } + } if (dump_) { - for (unsigned int i = 0; i < filteredMetaPrimitives.size(); i++) { - stringstream ss; - ss << iEvent.id().event() << " filtered mp " << i << ": "; - printmP(ss.str(), filteredMetaPrimitives.at(i)); + for (auto& ch_filteredMetaPrimitives : filteredMetaPrimitives) { + for (unsigned int i = 0; i < ch_filteredMetaPrimitives.second.size(); i++) { + stringstream ss; + ss << iEvent.id().event() << " filtered mp " << i << ": "; + printmP(ss.str(), ch_filteredMetaPrimitives.second.at(i)); + } } } - metaPrimitives.clear(); - metaPrimitives.erase(metaPrimitives.begin(), metaPrimitives.end()); + skip_processing_ = skip_processing_ || output_slfilter_; + confirmedMetaPrimitives.clear(); if (debug_) - LogDebug("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " - << filteredMetaPrimitives.size() << " filteredMetaPrimitives (superlayer)" - << std::endl; + for (auto& ch_filteredMetaPrimitives : filteredMetaPrimitives) { + LogDebug("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " + << ch_filteredMetaPrimitives.second.size() << " filteredMetaPrimitives (superlayer)" + << std::endl; + } if (debug_) LogDebug("DTTrigPhase2Prod") << "filteredMetaPrimitives: starting correlations" << std::endl; @@ -443,64 +718,106 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { //// CORRELATION: ///////////////////////////////////// - std::vector correlatedMetaPrimitives; - if (algo_ == Standard) - mpathassociator_->run(iEvent, iEventSetup, dtdigis, filteredMetaPrimitives, correlatedMetaPrimitives); - else { - for (const auto& muonpath : outmpaths) { - correlatedMetaPrimitives.emplace_back(muonpath->rawId(), - (double)muonpath->bxTimeValue(), - muonpath->horizPos(), - muonpath->tanPhi(), - muonpath->phi(), - muonpath->phiB(), - muonpath->phi_cmssw(), - muonpath->phiB_cmssw(), - muonpath->chiSquare(), - (int)muonpath->quality(), - muonpath->primitive(0)->channelId(), - muonpath->primitive(0)->tdcTimeStamp(), - muonpath->primitive(0)->laterality(), - muonpath->primitive(1)->channelId(), - muonpath->primitive(1)->tdcTimeStamp(), - muonpath->primitive(1)->laterality(), - muonpath->primitive(2)->channelId(), - muonpath->primitive(2)->tdcTimeStamp(), - muonpath->primitive(2)->laterality(), - muonpath->primitive(3)->channelId(), - muonpath->primitive(3)->tdcTimeStamp(), - muonpath->primitive(3)->laterality(), - muonpath->primitive(4)->channelId(), - muonpath->primitive(4)->tdcTimeStamp(), - muonpath->primitive(4)->laterality(), - muonpath->primitive(5)->channelId(), - muonpath->primitive(5)->tdcTimeStamp(), - muonpath->primitive(5)->laterality(), - muonpath->primitive(6)->channelId(), - muonpath->primitive(6)->tdcTimeStamp(), - muonpath->primitive(6)->laterality(), - muonpath->primitive(7)->channelId(), - muonpath->primitive(7)->tdcTimeStamp(), - muonpath->primitive(7)->laterality()); + std::map> correlatedMetaPrimitives; + if (algo_ == Standard) { + for (auto& ch_filteredMetaPrimitives : filteredMetaPrimitives) { + if (!skip_processing_) + mpathassociator_->run(iEvent, + iEventSetup, + ch_filteredMetaPrimitives.second, + correlatedMetaPrimitives[ch_filteredMetaPrimitives.first]); + else + for (auto& mp : ch_filteredMetaPrimitives.second) { + correlatedMetaPrimitives[ch_filteredMetaPrimitives.first].push_back(mp); + } + } + } else { + for (auto& ch_outmpaths : outmpaths) { + for (const auto& muonpath : ch_outmpaths.second) { + correlatedMetaPrimitives[ch_outmpaths.first].emplace_back(muonpath->rawId(), + (double)muonpath->bxTimeValue(), + muonpath->horizPos(), + muonpath->tanPhi(), + muonpath->phi(), + muonpath->phiB(), + muonpath->phi_cmssw(), + muonpath->phiB_cmssw(), + muonpath->chiSquare(), + (int)muonpath->quality(), + muonpath->primitive(0)->channelId(), + muonpath->primitive(0)->tdcTimeStamp(), + muonpath->primitive(0)->laterality(), + muonpath->primitive(1)->channelId(), + muonpath->primitive(1)->tdcTimeStamp(), + muonpath->primitive(1)->laterality(), + muonpath->primitive(2)->channelId(), + muonpath->primitive(2)->tdcTimeStamp(), + muonpath->primitive(2)->laterality(), + muonpath->primitive(3)->channelId(), + muonpath->primitive(3)->tdcTimeStamp(), + muonpath->primitive(3)->laterality(), + muonpath->primitive(4)->channelId(), + muonpath->primitive(4)->tdcTimeStamp(), + muonpath->primitive(4)->laterality(), + muonpath->primitive(5)->channelId(), + muonpath->primitive(5)->tdcTimeStamp(), + muonpath->primitive(5)->laterality(), + muonpath->primitive(6)->channelId(), + muonpath->primitive(6)->tdcTimeStamp(), + muonpath->primitive(6)->laterality(), + muonpath->primitive(7)->channelId(), + muonpath->primitive(7)->tdcTimeStamp(), + muonpath->primitive(7)->laterality()); + } } } - filteredMetaPrimitives.clear(); - if (debug_) - LogDebug("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " - << correlatedMetaPrimitives.size() << " correlatedMetPrimitives (chamber)"; + skip_processing_ = skip_processing_ || output_matcher_; + if (debug_) + for (auto& ch_correlatedMetaPrimitives : correlatedMetaPrimitives) { + LogDebug("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " + << ch_correlatedMetaPrimitives.second.size() << " correlatedMetPrimitives (chamber)"; + } if (dump_) { - LogInfo("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " - << correlatedMetaPrimitives.size() << " correlatedMetPrimitives (chamber)"; + for (auto& ch_correlatedMetaPrimitives : correlatedMetaPrimitives) { + LogDebug("DTTrigPhase2Prod") << "DTp2 in event:" << iEvent.id().event() << " we found " + << ch_correlatedMetaPrimitives.second.size() << " correlatedMetPrimitives (chamber)"; + } + for (auto& ch_correlatedMetaPrimitives : correlatedMetaPrimitives) { + for (unsigned int i = 0; i < ch_correlatedMetaPrimitives.second.size(); i++) { + stringstream ss; + ss << iEvent.id().event() << " correlated mp " << i << ": "; + printmPC(ss.str(), ch_correlatedMetaPrimitives.second.at(i)); + } + } + } - for (unsigned int i = 0; i < correlatedMetaPrimitives.size(); i++) { - stringstream ss; - ss << iEvent.id().event() << " correlated mp " << i << ": "; - printmPC(ss.str(), correlatedMetaPrimitives.at(i)); + // Correlated Filtering + std::map> filtCorrelatedMetaPrimitives; + if (algo_ == Standard) { + for (auto& ch_filteredMetaPrimitives : filteredMetaPrimitives) { + if (!skip_processing_) + mpathcorfilter_->run(iEvent, + iEventSetup, + ch_filteredMetaPrimitives.second, + correlatedMetaPrimitives[ch_filteredMetaPrimitives.first], + filtCorrelatedMetaPrimitives[ch_filteredMetaPrimitives.first]); + else { + for (auto& mp : ch_filteredMetaPrimitives.second) { + filtCorrelatedMetaPrimitives[ch_filteredMetaPrimitives.first].push_back(mp); + } + if (output_matcher_) + for (auto& mp : correlatedMetaPrimitives[ch_filteredMetaPrimitives.first]) { + filtCorrelatedMetaPrimitives[ch_filteredMetaPrimitives.first].push_back(mp); + } + } } } + correlatedMetaPrimitives.clear(); + filteredMetaPrimitives.clear(); + double shift_back = 0; if (scenario_ == MC) //scope for MC shift_back = 400; @@ -513,7 +830,9 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { if (useRPC_) { rpc_integrator_->initialise(iEventSetup, shift_back); rpc_integrator_->prepareMetaPrimitives(rpcRecHits); - rpc_integrator_->matchWithDTAndUseRPCTime(correlatedMetaPrimitives); + for (auto& ch_correlatedMetaPrimitives : filtCorrelatedMetaPrimitives) { + rpc_integrator_->matchWithDTAndUseRPCTime(ch_correlatedMetaPrimitives.second); // Probably this is a FIXME + } rpc_integrator_->makeRPCOnlySegments(); rpc_integrator_->storeRPCSingleHits(); rpc_integrator_->removeRPCHitsUsed(); @@ -526,148 +845,157 @@ void DTTrigPhase2Prod::produce(Event& iEvent, const EventSetup& iEventSetup) { vector outExtP2Th; // Assigning index value - assignIndex(correlatedMetaPrimitives); - for (const auto& metaPrimitiveIt : correlatedMetaPrimitives) { - DTChamberId chId(metaPrimitiveIt.rawId); - DTSuperLayerId slId(metaPrimitiveIt.rawId); - if (debug_) - LogDebug("DTTrigPhase2Prod") << "looping in final vector: SuperLayerId" << chId << " x=" << metaPrimitiveIt.x - << " quality=" << metaPrimitiveIt.quality - << " BX=" << round(metaPrimitiveIt.t0 / 25.) << " index=" << metaPrimitiveIt.index; - - int sectorTP = chId.sector(); - //sectors 13 and 14 exist only for the outermost stations for sectors 4 and 10 respectively - //due to the larger MB4 that are divided into two. - if (sectorTP == 13) - sectorTP = 4; - if (sectorTP == 14) - sectorTP = 10; - sectorTP = sectorTP - 1; - int sl = 0; - if (metaPrimitiveIt.quality < LOWLOWQ || metaPrimitiveIt.quality == CHIGHQ) { - if (inner(metaPrimitiveIt)) - sl = 1; - else - sl = 3; + if (!skip_processing_) + for (auto& ch_correlatedMetaPrimitives : filtCorrelatedMetaPrimitives) { + assignIndex(ch_correlatedMetaPrimitives.second); } - if (debug_) - LogDebug("DTTrigPhase2Prod") << "pushing back phase-2 dataformat carlo-federica dataformat"; - - if (slId.superLayer() != 2) { - if (df_extended_ == 1 || df_extended_ == 2) { - int pathWireId[8] = {metaPrimitiveIt.wi1, - metaPrimitiveIt.wi2, - metaPrimitiveIt.wi3, - metaPrimitiveIt.wi4, - metaPrimitiveIt.wi5, - metaPrimitiveIt.wi6, - metaPrimitiveIt.wi7, - metaPrimitiveIt.wi8}; - - int pathTDC[8] = {max((int)round(metaPrimitiveIt.tdc1 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc2 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc3 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc4 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc5 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc6 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc7 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc8 - shift_back * LHC_CLK_FREQ), -1)}; - - int pathLat[8] = {metaPrimitiveIt.lat1, - metaPrimitiveIt.lat2, - metaPrimitiveIt.lat3, - metaPrimitiveIt.lat4, - metaPrimitiveIt.lat5, - metaPrimitiveIt.lat6, - metaPrimitiveIt.lat7, - metaPrimitiveIt.lat8}; - - // phiTP (extended DF) - outExtP2Ph.emplace_back( - L1Phase2MuDTExtPhDigi((int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, - chId.wheel(), // uwh (m_wheel) - sectorTP, // usc (m_sector) - chId.station(), // ust (m_station) - sl, // ust (m_station) - (int)round(metaPrimitiveIt.phi * PHIRES_CONV), // uphi (m_phiAngle) - (int)round(metaPrimitiveIt.phiB * PHIBRES_CONV), // uphib (m_phiBending) - metaPrimitiveIt.quality, // uqua (m_qualityCode) - metaPrimitiveIt.index, // uind (m_segmentIndex) - (int)round(metaPrimitiveIt.t0) - shift_back * LHC_CLK_FREQ, // ut0 (m_t0Segment) - (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) - (int)round(metaPrimitiveIt.x * 1000), // ux (m_xLocal) - (int)round(metaPrimitiveIt.tanPhi * 1000), // utan (m_tanPsi) - (int)round(metaPrimitiveIt.phi_cmssw * PHIRES_CONV), // uphi (m_phiAngleCMSSW) - (int)round(metaPrimitiveIt.phiB_cmssw * PHIBRES_CONV), // uphib (m_phiBendingCMSSW) - metaPrimitiveIt.rpcFlag, // urpc (m_rpcFlag) - pathWireId, - pathTDC, - pathLat)); - } - if (df_extended_ == 0 || df_extended_ == 2) { - // phiTP (standard DF) - outP2Ph.push_back(L1Phase2MuDTPhDigi( - (int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, - chId.wheel(), // uwh (m_wheel) - sectorTP, // usc (m_sector) - chId.station(), // ust (m_station) - sl, // ust (m_station) - (int)round(metaPrimitiveIt.phi * PHIRES_CONV), // uphi (_phiAngle) - (int)round(metaPrimitiveIt.phiB * PHIBRES_CONV), // uphib (m_phiBending) - metaPrimitiveIt.quality, // uqua (m_qualityCode) - metaPrimitiveIt.index, // uind (m_segmentIndex) - (int)round(metaPrimitiveIt.t0) - shift_back * LHC_CLK_FREQ, // ut0 (m_t0Segment) - (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) - metaPrimitiveIt.rpcFlag // urpc (m_rpcFlag) - )); + for (auto& ch_correlatedMetaPrimitives : filtCorrelatedMetaPrimitives) { + for (const auto& metaPrimitiveIt : ch_correlatedMetaPrimitives.second) { + DTChamberId chId(metaPrimitiveIt.rawId); + DTSuperLayerId slId(metaPrimitiveIt.rawId); + if (debug_) + LogDebug("DTTrigPhase2Prod") << "looping in final vector: SuperLayerId" << chId << " x=" << metaPrimitiveIt.x + << " quality=" << metaPrimitiveIt.quality + << " BX=" << round(metaPrimitiveIt.t0 / 25.) << " index=" << metaPrimitiveIt.index; + + int sectorTP = chId.sector(); + //sectors 13 and 14 exist only for the outermost stations for sectors 4 and 10 respectively + //due to the larger MB4 that are divided into two. + if (sectorTP == 13) + sectorTP = 4; + if (sectorTP == 14) + sectorTP = 10; + sectorTP = sectorTP - 1; + int sl = 0; + if (metaPrimitiveIt.quality < LOWLOWQ || metaPrimitiveIt.quality == CHIGHQ) { + if (inner(metaPrimitiveIt)) + sl = 1; + else + sl = 3; } - } else { - if (df_extended_ == 1 || df_extended_ == 2) { - int pathWireId[4] = {metaPrimitiveIt.wi1, metaPrimitiveIt.wi2, metaPrimitiveIt.wi3, metaPrimitiveIt.wi4}; - - int pathTDC[4] = {max((int)round(metaPrimitiveIt.tdc1 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc2 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc3 - shift_back * LHC_CLK_FREQ), -1), - max((int)round(metaPrimitiveIt.tdc4 - shift_back * LHC_CLK_FREQ), -1)}; - - int pathLat[4] = {metaPrimitiveIt.lat1, metaPrimitiveIt.lat2, metaPrimitiveIt.lat3, metaPrimitiveIt.lat4}; - - // thTP (extended DF) - outExtP2Th.emplace_back( - L1Phase2MuDTExtThDigi((int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, - chId.wheel(), // uwh (m_wheel) - sectorTP, // usc (m_sector) - chId.station(), // ust (m_station) - (int)round(metaPrimitiveIt.phi * ZRES_CONV), // uz (m_zGlobal) - (int)round(metaPrimitiveIt.phiB * KRES_CONV), // uk (m_kSlope) - metaPrimitiveIt.quality, // uqua (m_qualityCode) - metaPrimitiveIt.index, // uind (m_segmentIndex) - (int)round(metaPrimitiveIt.t0) - shift_back * LHC_CLK_FREQ, // ut0 (m_t0Segment) - (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) - (int)round(metaPrimitiveIt.x * 1000), // ux (m_yLocal) - (int)round(metaPrimitiveIt.phi_cmssw * PHIRES_CONV), // uphi (m_zCMSSW) - (int)round(metaPrimitiveIt.phiB_cmssw * PHIBRES_CONV), // uphib (m_kCMSSW) - metaPrimitiveIt.rpcFlag, // urpc (m_rpcFlag) - pathWireId, - pathTDC, - pathLat)); - } - if (df_extended_ == 0 || df_extended_ == 2) { - // thTP (standard DF) - outP2Th.push_back(L1Phase2MuDTThDigi( - (int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, - chId.wheel(), // uwh (m_wheel) - sectorTP, // usc (m_sector) - chId.station(), // ust (m_station) - (int)round(metaPrimitiveIt.phi * ZRES_CONV), // uz (m_zGlobal) - (int)round(metaPrimitiveIt.phiB * KRES_CONV), // uk (m_kSlope) - metaPrimitiveIt.quality, // uqua (m_qualityCode) - metaPrimitiveIt.index, // uind (m_segmentIndex) - (int)round(metaPrimitiveIt.t0) - shift_back * LHC_CLK_FREQ, // ut0 (m_t0Segment) - (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) - metaPrimitiveIt.rpcFlag // urpc (m_rpcFlag) - )); + + float tp_t0 = + (metaPrimitiveIt.t0 - shift_back * LHC_CLK_FREQ) * ((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ); + + if (debug_) + LogDebug("DTTrigPhase2Prod") << "pushing back phase-2 dataformat carlo-federica dataformat"; + + if (slId.superLayer() != 2) { + if (df_extended_ == 1 || df_extended_ == 2) { + int pathWireId[8] = {metaPrimitiveIt.wi1, + metaPrimitiveIt.wi2, + metaPrimitiveIt.wi3, + metaPrimitiveIt.wi4, + metaPrimitiveIt.wi5, + metaPrimitiveIt.wi6, + metaPrimitiveIt.wi7, + metaPrimitiveIt.wi8}; + + int pathTDC[8] = {max((int)round(metaPrimitiveIt.tdc1 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc2 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc3 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc4 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc5 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc6 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc7 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc8 - shift_back * LHC_CLK_FREQ), -1)}; + + int pathLat[8] = {metaPrimitiveIt.lat1, + metaPrimitiveIt.lat2, + metaPrimitiveIt.lat3, + metaPrimitiveIt.lat4, + metaPrimitiveIt.lat5, + metaPrimitiveIt.lat6, + metaPrimitiveIt.lat7, + metaPrimitiveIt.lat8}; + + // phiTP (extended DF) + outExtP2Ph.emplace_back( + L1Phase2MuDTExtPhDigi((int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, + chId.wheel(), // uwh (m_wheel) + sectorTP, // usc (m_sector) + chId.station(), // ust (m_station) + sl, // ust (m_station) + (int)round(metaPrimitiveIt.phi * PHIRES_CONV), // uphi (m_phiAngle) + (int)round(metaPrimitiveIt.phiB * PHIBRES_CONV), // uphib (m_phiBending) + metaPrimitiveIt.quality, // uqua (m_qualityCode) + metaPrimitiveIt.index, // uind (m_segmentIndex) + tp_t0, // ut0 (m_t0Segment) + (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) + (int)round(metaPrimitiveIt.x * 1000), // ux (m_xLocal) + (int)round(metaPrimitiveIt.tanPhi * 1000), // utan (m_tanPsi) + (int)round(metaPrimitiveIt.phi_cmssw * PHIRES_CONV), // uphi (m_phiAngleCMSSW) + (int)round(metaPrimitiveIt.phiB_cmssw * PHIBRES_CONV), // uphib (m_phiBendingCMSSW) + metaPrimitiveIt.rpcFlag, // urpc (m_rpcFlag) + pathWireId, + pathTDC, + pathLat)); + } + if (df_extended_ == 0 || df_extended_ == 2) { + // phiTP (standard DF) + outP2Ph.push_back(L1Phase2MuDTPhDigi( + (int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, + chId.wheel(), // uwh (m_wheel) + sectorTP, // usc (m_sector) + chId.station(), // ust (m_station) + sl, // ust (m_station) + (int)round(metaPrimitiveIt.phi * PHIRES_CONV), // uphi (_phiAngle) + (int)round(metaPrimitiveIt.phiB * PHIBRES_CONV), // uphib (m_phiBending) + metaPrimitiveIt.quality, // uqua (m_qualityCode) + metaPrimitiveIt.index, // uind (m_segmentIndex) + tp_t0, // ut0 (m_t0Segment) + (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) + metaPrimitiveIt.rpcFlag // urpc (m_rpcFlag) + )); + } + } else { + if (df_extended_ == 1 || df_extended_ == 2) { + int pathWireId[4] = {metaPrimitiveIt.wi1, metaPrimitiveIt.wi2, metaPrimitiveIt.wi3, metaPrimitiveIt.wi4}; + + int pathTDC[4] = {max((int)round(metaPrimitiveIt.tdc1 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc2 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc3 - shift_back * LHC_CLK_FREQ), -1), + max((int)round(metaPrimitiveIt.tdc4 - shift_back * LHC_CLK_FREQ), -1)}; + + int pathLat[4] = {metaPrimitiveIt.lat1, metaPrimitiveIt.lat2, metaPrimitiveIt.lat3, metaPrimitiveIt.lat4}; + + // thTP (extended DF) + outExtP2Th.emplace_back( + L1Phase2MuDTExtThDigi((int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, + chId.wheel(), // uwh (m_wheel) + sectorTP, // usc (m_sector) + chId.station(), // ust (m_station) + (int)round(metaPrimitiveIt.phi * ZRES_CONV), // uz (m_zGlobal) + (int)round(metaPrimitiveIt.phiB * KRES_CONV), // uk (m_kSlope) + metaPrimitiveIt.quality, // uqua (m_qualityCode) + metaPrimitiveIt.index, // uind (m_segmentIndex) + tp_t0, // ut0 (m_t0Segment) + (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) + (int)round(metaPrimitiveIt.x * 1000), // ux (m_yLocal) + (int)round(metaPrimitiveIt.phi_cmssw * ZRES_CONV), // uphi (m_zCMSSW) + (int)round(metaPrimitiveIt.phiB_cmssw * KRES_CONV), // uphib (m_kCMSSW) + metaPrimitiveIt.rpcFlag, // urpc (m_rpcFlag) + pathWireId, + pathTDC, + pathLat)); + } + if (df_extended_ == 0 || df_extended_ == 2) { + // thTP (standard DF) + outP2Th.push_back(L1Phase2MuDTThDigi( + (int)round(metaPrimitiveIt.t0 / (float)LHC_CLK_FREQ) - shift_back, + chId.wheel(), // uwh (m_wheel) + sectorTP, // usc (m_sector) + chId.station(), // ust (m_station) + (int)round(metaPrimitiveIt.phi * ZRES_CONV), // uz (m_zGlobal) + (int)round(metaPrimitiveIt.phiB * KRES_CONV), // uk (m_kSlope) + metaPrimitiveIt.quality, // uqua (m_qualityCode) + metaPrimitiveIt.index, // uind (m_segmentIndex) + tp_t0, // ut0 (m_t0Segment) + (int)round(metaPrimitiveIt.chi2 * CHI2RES_CONV), // uchi2 (m_chi2Segment) + metaPrimitiveIt.rpcFlag // urpc (m_rpcFlag) + )); + } } } } @@ -747,6 +1075,16 @@ void DTTrigPhase2Prod::printmP(const string& ss, const metaPrimitive& mP) const << setw(13) << left << mP.chi2 << " r:" << rango(mP); } +void DTTrigPhase2Prod::printmP(const metaPrimitive& mP) const { + DTSuperLayerId slId(mP.rawId); + LogInfo("DTTrigPhase2Prod") << (int)slId << "\t " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 + << " " << setw(2) << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(5) + << left << mP.tdc1 << " " << setw(5) << left << mP.tdc2 << " " << setw(5) << left + << mP.tdc3 << " " << setw(5) << left << mP.tdc4 << " " << setw(10) << right << mP.x << " " + << setw(9) << left << mP.tanPhi << " " << setw(5) << left << mP.t0 << " " << setw(13) + << left << mP.chi2 << " r:" << rango(mP) << std::endl; +} + void DTTrigPhase2Prod::printmPC(const string& ss, const metaPrimitive& mP) const { DTChamberId ChId(mP.rawId); LogInfo("DTTrigPhase2Prod") << ss << (int)ChId << "\t " << setw(2) << left << mP.wi1 << " " << setw(2) << left @@ -764,6 +1102,23 @@ void DTTrigPhase2Prod::printmPC(const string& ss, const metaPrimitive& mP) const << left << mP.chi2 << " r:" << rango(mP); } +void DTTrigPhase2Prod::printmPC(const metaPrimitive& mP) const { + DTChamberId ChId(mP.rawId); + LogInfo("DTTrigPhase2Prod") << (int)ChId << "\t " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 + << " " << setw(2) << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(2) + << left << mP.wi5 << " " << setw(2) << left << mP.wi6 << " " << setw(2) << left << mP.wi7 + << " " << setw(2) << left << mP.wi8 << " " << setw(5) << left << mP.tdc1 << " " << setw(5) + << left << mP.tdc2 << " " << setw(5) << left << mP.tdc3 << " " << setw(5) << left + << mP.tdc4 << " " << setw(5) << left << mP.tdc5 << " " << setw(5) << left << mP.tdc6 + << " " << setw(5) << left << mP.tdc7 << " " << setw(5) << left << mP.tdc8 << " " + << setw(2) << left << mP.lat1 << " " << setw(2) << left << mP.lat2 << " " << setw(2) + << left << mP.lat3 << " " << setw(2) << left << mP.lat4 << " " << setw(2) << left + << mP.lat5 << " " << setw(2) << left << mP.lat6 << " " << setw(2) << left << mP.lat7 + << " " << setw(2) << left << mP.lat8 << " " << setw(10) << right << mP.x << " " << setw(9) + << left << mP.tanPhi << " " << setw(5) << left << mP.t0 << " " << setw(13) << left + << mP.chi2 << " r:" << rango(mP) << std::endl; +} + int DTTrigPhase2Prod::rango(const metaPrimitive& mp) const { if (mp.quality == 1 or mp.quality == 2) return 3; @@ -880,7 +1235,6 @@ void DTTrigPhase2Prod::fillDescriptions(edm::ConfigurationDescriptions& descript // dtTriggerPhase2PrimitiveDigis edm::ParameterSetDescription desc; desc.add("digiTag", edm::InputTag("CalibratedDigis")); - desc.add("trigger_with_sl", 4); desc.add("timeTolerance", 999999); desc.add("tanPhiTh", 1.0); desc.add("tanPhiThw2max", 1.3); @@ -900,12 +1254,26 @@ void DTTrigPhase2Prod::fillDescriptions(edm::ConfigurationDescriptions& descript desc.add("scenario", 0); desc.add("df_extended", 0); desc.add("max_primitives", 999); + desc.add("output_mixer", false); + desc.add("output_latpredictor", false); + desc.add("output_slfitter", false); + desc.add("output_slfilter", false); + desc.add("output_confirmed", false); + desc.add("output_matcher", false); desc.add("ttrig_filename", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/wire_rawId_ttrig.txt")); desc.add("z_filename", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/wire_rawId_z.txt")); + desc.add("lut_sl1", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/fitterlut_sl1.dat")); + desc.add("lut_sl2", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/fitterlut_slx.dat")); + desc.add("lut_sl3", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/fitterlut_sl3.dat")); + desc.add("lut_2sl", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/fitterlut_2sl.dat")); desc.add("shift_filename", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/wire_rawId_x.txt")); + desc.add("maxdrift_filename", + edm::FileInPath("L1Trigger/DTTriggerPhase2/data/drift_time_per_chamber.txt")); desc.add("shift_theta_filename", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/theta_shift.txt")); desc.add("global_coords_filename", edm::FileInPath("L1Trigger/DTTriggerPhase2/data/global_coord_perp_x_phi0.txt")); + desc.add("laterality_filename", + edm::FileInPath("L1Trigger/DTTriggerPhase2/data/lat_predictions.dat")); desc.add("algo", 0); desc.add("minHits4Fit", 3); desc.add("splitPathPerSL", true); diff --git a/L1Trigger/DTTriggerPhase2/python/dtTriggerPhase2PrimitiveDigis_cfi.py b/L1Trigger/DTTriggerPhase2/python/dtTriggerPhase2PrimitiveDigis_cfi.py index 5e3f920736e03..d044ecf89923b 100644 --- a/L1Trigger/DTTriggerPhase2/python/dtTriggerPhase2PrimitiveDigis_cfi.py +++ b/L1Trigger/DTTriggerPhase2/python/dtTriggerPhase2PrimitiveDigis_cfi.py @@ -6,7 +6,6 @@ dtTriggerPhase2PrimitiveDigis = cms.EDProducer("DTTrigPhase2Prod", digiTag = cms.InputTag("CalibratedDigis"), - trigger_with_sl = cms.int32(4), tanPhiTh = cms.double(1.), tanPhiThw2max = cms.double(1.3), tanPhiThw2min = cms.double(0.5), @@ -26,11 +25,24 @@ df_extended = cms.int32(0), # DF: 0 for standard, 1 for extended, 2 for both max_primitives = cms.int32(999), + output_mixer = cms.bool(False), + output_latpredictor = cms.bool(False), + output_slfitter = cms.bool(False), + output_slfilter = cms.bool(False), + output_confirmed = cms.bool(False), + output_matcher = cms.bool(False), + ttrig_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/wire_rawId_ttrig.txt'), z_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/wire_rawId_z.txt'), + lut_sl1 = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/fitterlut_sl1.dat'), + lut_sl2 = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/fitterlut_slx.dat'), + lut_sl3 = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/fitterlut_sl3.dat'), + lut_2sl = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/fitterlut_2sl.dat'), shift_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/wire_rawId_x.txt'), shift_theta_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/theta_shift.txt'), + maxdrift_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/drift_time_per_chamber.txt'), global_coords_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/global_coord_perp_x_phi0.txt'), + laterality_filename = cms.FileInPath('L1Trigger/DTTriggerPhase2/data/lat_predictions.dat'), algo = cms.int32(0), # 0 = STD gr., 2 = Hough transform, 1 = PseudoBayes Approach minHits4Fit = cms.int32(3), diff --git a/L1Trigger/DTTriggerPhase2/src/LateralityBasicProvider.cc b/L1Trigger/DTTriggerPhase2/src/LateralityBasicProvider.cc new file mode 100644 index 0000000000000..a46a2e82b5e87 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/LateralityBasicProvider.cc @@ -0,0 +1,100 @@ +#include "L1Trigger/DTTriggerPhase2/interface/LateralityBasicProvider.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; +// ============================================================================ +// Constructors and destructor +// ============================================================================ +LateralityBasicProvider::LateralityBasicProvider(const ParameterSet &pset, edm::ConsumesCollector &iC) + : LateralityProvider(pset, iC), debug_(pset.getUntrackedParameter("debug")) { + if (debug_) + LogDebug("LateralityBasicProvider") << "LateralityBasicProvider: constructor"; + + fill_lat_combinations(); +} + +LateralityBasicProvider::~LateralityBasicProvider() { + if (debug_) + LogDebug("LateralityBasicProvider") << "LateralityBasicProvider: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void LateralityBasicProvider::initialise(const edm::EventSetup &iEventSetup) { + if (debug_) + LogDebug("LateralityBasicProvider") << "LateralityBasicProvider::initialiase"; +} + +void LateralityBasicProvider::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &muonpaths, + std::vector &lateralities) { + if (debug_) + LogDebug("LateralityBasicProvider") << "LateralityBasicProvider: run"; + + // fit per SL (need to allow for multiple outputs for a single mpath) + for (auto &muonpath : muonpaths) { + analyze(muonpath, lateralities); + } +} + +void LateralityBasicProvider::finish() { + if (debug_) + LogDebug("LateralityBasicProvider") << "LateralityBasicProvider: finish"; +}; + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +void LateralityBasicProvider::analyze(MuonPathPtr &inMPath, std::vector &lateralities) { + if (debug_) + LogDebug("LateralityBasicProvider") << "DTp2:analyze \t\t\t\t starts"; + for (auto &lat_combination : lat_combinations) { + if (inMPath->missingLayer() == lat_combination.missing_layer && + inMPath->cellLayout()[0] == lat_combination.cellLayout[0] && + inMPath->cellLayout()[1] == lat_combination.cellLayout[1] && + inMPath->cellLayout()[2] == lat_combination.cellLayout[2] && + inMPath->cellLayout()[3] == lat_combination.cellLayout[3]) { + lateralities.push_back(lat_combination.latcombs); + return; + } + } + lateralities.push_back(LAT_VECTOR_NULL); + return; +} + +void LateralityBasicProvider::fill_lat_combinations() { + lat_combinations.push_back({-1, {0, 0, 0, -1}, {{0, 0, 0, 1}, {0, 0, 1, 1}, {0, 1, 1, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({-1, {0, 0, 1, -1}, {{0, 0, 1, 0}, {0, 1, 1, 0}, {1, 1, 1, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({-1, {0, 1, 0, -1}, {{0, 1, 0, 0}, {0, 1, 0, 1}, {1, 1, 0, 0}, {1, 1, 0, 1}}}); + lat_combinations.push_back({-1, {0, 1, 1, -1}, {{0, 1, 0, 0}, {0, 1, 1, 0}, {0, 1, 1, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({-1, {1, 0, 0, -1}, {{1, 0, 0, 0}, {1, 0, 0, 1}, {1, 0, 1, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({-1, {1, 0, 1, -1}, {{0, 0, 1, 0}, {0, 0, 1, 1}, {1, 0, 1, 0}, {1, 0, 1, 1}}}); + lat_combinations.push_back({-1, {1, 1, 0, -1}, {{0, 0, 0, 1}, {1, 0, 0, 1}, {1, 1, 0, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({-1, {1, 1, 1, -1}, {{1, 0, 0, 0}, {1, 1, 0, 0}, {1, 1, 1, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({0, {0, 0, 0, -1}, {{0, 0, 0, 1}, {0, 0, 1, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({0, {0, 0, 1, -1}, {{0, 0, 1, 0}, {0, 0, 1, 1}, {0, 1, 1, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({0, {0, 1, 0, -1}, {{0, 0, 0, 1}, {0, 1, 0, 0}, {0, 1, 0, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({0, {0, 1, 1, -1}, {{0, 1, 0, 0}, {0, 1, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {0, 0, 0, -1}, {{0, 0, 0, 1}, {0, 0, 1, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {0, 0, 1, -1}, {{0, 0, 1, 0}, {1, 0, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {0, 1, 0, -1}, {{0, 0, 0, 1}, {1, 0, 0, 0}, {1, 0, 0, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {0, 1, 1, -1}, {{0, 0, 1, 0}, {0, 0, 1, 1}, {1, 0, 1, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {1, 1, 0, -1}, {{0, 0, 0, 1}, {1, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({1, {1, 1, 1, -1}, {{1, 0, 0, 0}, {1, 0, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {0, 0, 0, -1}, {{0, 0, 0, 1}, {0, 1, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {0, 0, 1, -1}, {{0, 1, 0, 0}, {0, 1, 0, 1}, {1, 1, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {0, 1, 1, -1}, {{0, 1, 0, 0}, {0, 1, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {1, 0, 0, -1}, {{1, 0, 0, 0}, {1, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {1, 0, 1, -1}, {{0, 0, 0, 1}, {1, 0, 0, 0}, {1, 0, 0, 1}, {0, 0, 0, 0}}}); + lat_combinations.push_back({2, {1, 1, 1, -1}, {{1, 0, 0, 0}, {1, 1, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({3, {0, 0, 0, -1}, {{0, 0, 1, 0}, {0, 1, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({3, {0, 1, 0, -1}, {{0, 1, 0, 0}, {0, 1, 1, 0}, {1, 1, 0, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({3, {1, 0, 0, -1}, {{0, 0, 1, 0}, {1, 0, 0, 0}, {1, 0, 1, 0}, {0, 0, 0, 0}}}); + lat_combinations.push_back({3, {1, 1, 0, -1}, {{1, 0, 0, 0}, {1, 1, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}}); +}; diff --git a/L1Trigger/DTTriggerPhase2/src/LateralityCoarsedProvider.cc b/L1Trigger/DTTriggerPhase2/src/LateralityCoarsedProvider.cc new file mode 100644 index 0000000000000..e37281d999e24 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/LateralityCoarsedProvider.cc @@ -0,0 +1,259 @@ +#include "L1Trigger/DTTriggerPhase2/interface/LateralityCoarsedProvider.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; +// ============================================================================ +// Constructors and destructor +// ============================================================================ +LateralityCoarsedProvider::LateralityCoarsedProvider(const ParameterSet &pset, edm::ConsumesCollector &iC) + : LateralityProvider(pset, iC), + debug_(pset.getUntrackedParameter("debug")), + laterality_filename_(pset.getParameter("laterality_filename")) { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "LateralityCoarsedProvider: constructor"; + + fill_lat_combinations(); +} + +LateralityCoarsedProvider::~LateralityCoarsedProvider() { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "LateralityCoarsedProvider: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void LateralityCoarsedProvider::initialise(const edm::EventSetup &iEventSetup) { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "LateralityCoarsedProvider::initialiase"; +} + +void LateralityCoarsedProvider::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &muonpaths, + std::vector &lateralities) { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "LateralityCoarsedProvider: run"; + + // fit per SL (need to allow for multiple outputs for a single mpath) + for (auto &muonpath : muonpaths) { + analyze(muonpath, lateralities); + } +} + +void LateralityCoarsedProvider::finish() { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "LateralityCoarsedProvider: finish"; +}; + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +void LateralityCoarsedProvider::analyze(MuonPathPtr &inMPath, std::vector &lateralities) { + if (debug_) + LogDebug("LateralityCoarsedProvider") << "DTp2:analyze \t\t\t\t starts"; + + auto coarsified_times = coarsify_times(inMPath); + + for (auto &lat_combination : lat_combinations) { + if (inMPath->missingLayer() == lat_combination.missing_layer && + inMPath->cellLayout()[0] == lat_combination.cellLayout[0] && + inMPath->cellLayout()[1] == lat_combination.cellLayout[1] && + inMPath->cellLayout()[2] == lat_combination.cellLayout[2] && + inMPath->cellLayout()[3] == lat_combination.cellLayout[3] && + coarsified_times[0] == lat_combination.coarsed_times[0] && + coarsified_times[1] == lat_combination.coarsed_times[1] && + coarsified_times[2] == lat_combination.coarsed_times[2] && + coarsified_times[3] == lat_combination.coarsed_times[3]) { + lateralities.push_back(lat_combination.latcombs); + return; + } + } + lateralities.push_back(LAT_VECTOR_NULL); + return; +} + +std::vector LateralityCoarsedProvider::coarsify_times(MuonPathPtr &inMPath) { + int max_time = -999; + // obtain the maximum time to do the coarsification + for (int layer = 0; layer < cmsdt::NUM_LAYERS; layer++) { + if (inMPath->missingLayer() == layer) + continue; + if (inMPath->primitive(layer)->tdcTimeStamp() > max_time) + max_time = inMPath->primitive(layer)->tdcTimeStamp(); + } + + // do the coarsification + std::vector coarsified_times; + for (int layer = 0; layer < cmsdt::NUM_LAYERS; layer++) { + if (inMPath->missingLayer() == layer) { + coarsified_times.push_back(-1); + continue; + } + auto coarsified_time = max_time - inMPath->primitive(layer)->tdcTimeStamp(); + // transform into tdc counts + coarsified_time = (int)round(((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ) * coarsified_time); + // keep the LAT_MSB_BITS + coarsified_time = coarsified_time >> (LAT_TOTAL_BITS - LAT_MSB_BITS); + + if (inMPath->missingLayer() == -1) { // 4-hit candidates + if (coarsified_time <= LAT_P0_4H) + coarsified_times.push_back(0); + else if (coarsified_time <= LAT_P1_4H) + coarsified_times.push_back(1); + else if (coarsified_time <= LAT_P2_4H) + coarsified_times.push_back(2); + else + coarsified_times.push_back(3); + } else { // 3-hit candidates + if (coarsified_time <= LAT_P0_3H) + coarsified_times.push_back(0); + else if (coarsified_time <= LAT_P1_3H) + coarsified_times.push_back(1); + else if (coarsified_time <= LAT_P2_3H) + coarsified_times.push_back(2); + else + coarsified_times.push_back(3); + } + } + return coarsified_times; +} + +void LateralityCoarsedProvider::fill_lat_combinations() { + std::ifstream latFile(laterality_filename_.fullPath()); // Open file + if (latFile.fail()) { + throw cms::Exception("Missing Input File") + << "LateralityCoarsedProvider::fill_lat_combinations() - Cannot find " << laterality_filename_.fullPath(); + return; + } + + std::string line; + + short line_counter = 0; // Line counter + + // Bit masks for every parameter + int _12bitMask = 0xFFF; // 12 bits + int _layoutMask = 0xE00; // 3 bits + int _is4HitMask = 0x100; // 1 bit + int _coarsedMask = 0xFF; // 8 bits + int _layerMask = 0xC0; // 2 bits + + while (std::getline(latFile, line)) { + if (line == "000000000000") { + line_counter++; + continue; + } //skip zeros + + if (line.size() == 12) { + std::vector> transformedVector = convertString(line); + latcomb lat0 = { + transformedVector[0][0], transformedVector[0][1], transformedVector[0][2], transformedVector[0][3]}; + latcomb lat1 = { + transformedVector[1][0], transformedVector[1][1], transformedVector[1][2], transformedVector[1][3]}; + latcomb lat2 = { + transformedVector[2][0], transformedVector[2][1], transformedVector[2][2], transformedVector[2][3]}; + + //Transforming line number to binary + short address = line_counter & _12bitMask; // 12 bits + + short layout = + (address & _layoutMask) >> 9; //Doing AND and displacing 9 bits to the right to obtain 3 bits of layout + short is4Hit = (address & _is4HitMask) >> 8; + short coarsed = address & _coarsedMask; + + short bit1Layout = (layout & (1)); + short bit2Layout = (layout & (1 << 1)) >> 1; + short bit3Layout = (layout & (1 << 2)) >> 2; + + //Logic implementation + short missingLayer = -1; + short layout_comb[NUM_LAYERS] = {bit3Layout, bit2Layout, bit1Layout, -1}; + short coarsedTimes[NUM_LAYERS] = {0, 0, 0, 0}; + + if (is4Hit != 1) { //3 hit case + missingLayer = + (coarsed & _layerMask) >> 6; //Missing layer is given by the two most significative bits of coarsed vector + coarsedTimes[missingLayer] = -1; //Missing layer set to -1 + } + + // Filling coarsedTimes vector without the missing layer + if (missingLayer != -1) { + switch (missingLayer) { + case 0: + coarsedTimes[1] = (coarsed & 0x30) >> 4; + coarsedTimes[2] = (coarsed & 0x0C) >> 2; + coarsedTimes[3] = coarsed & 0x03; + lat0 = {-1, transformedVector[0][1], transformedVector[0][2], transformedVector[0][3]}; + lat1 = {-1, transformedVector[1][1], transformedVector[1][2], transformedVector[1][3]}; + lat2 = {-1, transformedVector[2][1], transformedVector[2][2], transformedVector[2][3]}; + break; + case 1: + coarsedTimes[0] = (coarsed & 0x30) >> 4; + coarsedTimes[2] = (coarsed & 0x0C) >> 2; + coarsedTimes[3] = coarsed & 0x03; + lat0 = {transformedVector[0][0], -1, transformedVector[0][2], transformedVector[0][3]}; + lat1 = {transformedVector[1][0], -1, transformedVector[1][2], transformedVector[1][3]}; + lat2 = {transformedVector[2][0], -1, transformedVector[2][2], transformedVector[2][3]}; + break; + case 2: + coarsedTimes[0] = (coarsed & 0x30) >> 4; + coarsedTimes[1] = (coarsed & 0x0C) >> 2; + coarsedTimes[3] = coarsed & 0x03; + lat0 = {transformedVector[0][0], transformedVector[0][1], -1, transformedVector[0][3]}; + lat1 = {transformedVector[1][0], transformedVector[1][1], -1, transformedVector[1][3]}; + lat2 = {transformedVector[2][0], transformedVector[2][1], -1, transformedVector[2][3]}; + break; + case 3: + coarsedTimes[0] = (coarsed & 0x30) >> 4; + coarsedTimes[1] = (coarsed & 0x0C) >> 2; + coarsedTimes[2] = coarsed & 0x03; + lat0 = {transformedVector[0][0], transformedVector[0][1], transformedVector[0][2], -1}; + lat1 = {transformedVector[1][0], transformedVector[1][1], transformedVector[1][2], -1}; + lat2 = {transformedVector[2][0], transformedVector[2][1], transformedVector[2][2], -1}; + break; + + default: + break; + } + + } else { //4 hit case + coarsedTimes[0] = (coarsed & 0xC0) >> 6; + coarsedTimes[1] = (coarsed & 0x30) >> 4; + coarsedTimes[2] = (coarsed & 0x0C) >> 2; + coarsedTimes[3] = coarsed & 0x03; + } + + lat_coarsed_combination lat_temp = {missingLayer, + {layout_comb[0], layout_comb[1], layout_comb[2], layout_comb[3]}, + {coarsedTimes[0], coarsedTimes[1], coarsedTimes[2], coarsedTimes[3]}, + {lat0, lat1, lat2}}; + lat_combinations.push_back(lat_temp); + + } else { //size different from 12 + std::cerr << "Error: line " << line_counter << " does not contain 12 bits." << std::endl; + } + line_counter++; + }; + + //closing lateralities file + latFile.close(); +}; + +// Function to convert a 12 bit string in a a vector of 4 bit vectors +std::vector> LateralityCoarsedProvider::convertString(std::string chain) { + std::vector> result; + + for (size_t i = 0; i < chain.size(); i += 4) { + std::vector group; + for (size_t j = 0; j < 4; j++) { + group.push_back(chain[i + j] - '0'); // Convert the character to integer + } + result.push_back(group); + } + + return result; +} diff --git a/L1Trigger/DTTriggerPhase2/src/LateralityProvider.cc b/L1Trigger/DTTriggerPhase2/src/LateralityProvider.cc new file mode 100644 index 0000000000000..ae5e8486c9242 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/LateralityProvider.cc @@ -0,0 +1,24 @@ +#include "L1Trigger/DTTriggerPhase2/interface/LateralityProvider.h" + +using namespace edm; +using namespace std; + +// ============================================================================ +// Constructors and destructor +// ============================================================================ +LateralityProvider::LateralityProvider(const ParameterSet& pset, edm::ConsumesCollector& iC) + : debug_(pset.getUntrackedParameter("debug")) {} + +LateralityProvider::~LateralityProvider() {} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void LateralityProvider::initialise(const edm::EventSetup& iEventSetup) {} + +void LateralityProvider::finish(){}; + +void LateralityProvider::run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + MuonPathPtrs& inMpath, + std::vector& lateralities){}; diff --git a/L1Trigger/DTTriggerPhase2/src/MPCorFilter.cc b/L1Trigger/DTTriggerPhase2/src/MPCorFilter.cc new file mode 100644 index 0000000000000..0feec5b478522 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MPCorFilter.cc @@ -0,0 +1,242 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MPCorFilter.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +using namespace edm; +using namespace std; +using namespace cmsdt; + +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MPCorFilter::MPCorFilter(const ParameterSet &pset) + : MPFilter(pset), debug_(pset.getUntrackedParameter("debug")) {} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void MPCorFilter::initialise(const edm::EventSetup &iEventSetup) {} + +void MPCorFilter::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inSLMPaths, + std::vector &inCorMPaths, + std::vector &outMPaths) { + if (debug_) + LogDebug("MPCorFilter") << "MPCorFilter: run"; + + std::vector SL1metaPrimitives; + std::vector SL2metaPrimitives; + std::vector SL3metaPrimitives; + std::vector CormetaPrimitives; + uint32_t sl1Id_rawid = -1, sl2Id_rawid = -1, sl3Id_rawid = -1; + if (!inSLMPaths.empty()) { + int dum_sl_rawid = inSLMPaths[0].rawId; + DTSuperLayerId dumSlId(dum_sl_rawid); + + max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1]; + DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector()); + DTSuperLayerId sl1Id(ChId.rawId(), 1); + sl1Id_rawid = sl1Id.rawId(); + DTSuperLayerId sl2Id(ChId.rawId(), 2); + sl2Id_rawid = sl2Id.rawId(); + DTSuperLayerId sl3Id(ChId.rawId(), 3); + sl3Id_rawid = sl3Id.rawId(); + + for (const auto &metaprimitiveIt : inSLMPaths) { + if (metaprimitiveIt.rawId == sl1Id_rawid) { + SL1metaPrimitives.push_back(metaprimitiveIt); + } else if (metaprimitiveIt.rawId == sl3Id_rawid) + SL3metaPrimitives.push_back(metaprimitiveIt); + else if (metaprimitiveIt.rawId == sl2Id_rawid) + SL2metaPrimitives.push_back(metaprimitiveIt); + } + } + auto filteredMPs = filter(SL1metaPrimitives, SL2metaPrimitives, SL3metaPrimitives, inCorMPaths); + for (auto &mp : filteredMPs) + outMPaths.push_back(mp); +} + +void MPCorFilter::finish(){}; + +/////////////////////////// +/// OTHER METHODS + +std::vector MPCorFilter::filter(std::vector SL1mps, + std::vector SL2mps, + std::vector SL3mps, + std::vector Cormps) { + std::map mp_valid_per_bx; + std::map imp_per_bx_sl1; + for (auto &mp : SL1mps) { + int BX = mp.t0 / 25; + if (mp_valid_per_bx.find(BX) == mp_valid_per_bx.end()) { + mp_valid_per_bx[BX] = valid_cor_tp_arr_t(12); + } + + if (imp_per_bx_sl1.find(BX) == imp_per_bx_sl1.end()) { + imp_per_bx_sl1[BX] = 0; + } + + auto coarsed = coarsify(mp, 1); + mp_valid_per_bx[BX][imp_per_bx_sl1[BX]] = valid_cor_tp_t({true, mp, coarsed[3], coarsed[4], coarsed[5]}); + imp_per_bx_sl1[BX] += 2; + } + std::map imp_per_bx_sl3; + for (auto &mp : SL3mps) { + int BX = mp.t0 / 25; + if (mp_valid_per_bx.find(BX) == mp_valid_per_bx.end()) { + mp_valid_per_bx[BX] = valid_cor_tp_arr_t(12); + } + + if (imp_per_bx_sl3.find(BX) == imp_per_bx_sl3.end()) { + imp_per_bx_sl3[BX] = 1; + } + + auto coarsed = coarsify(mp, 3); + mp_valid_per_bx[BX][imp_per_bx_sl3[BX]] = valid_cor_tp_t({true, mp, coarsed[3], coarsed[4], coarsed[5]}); + imp_per_bx_sl3[BX] += 2; + } + + for (auto &mp : Cormps) { + int BX = mp.t0 / 25; + if (mp_valid_per_bx.find(BX) == mp_valid_per_bx.end()) { + mp_valid_per_bx[BX] = valid_cor_tp_arr_t(12); + } + auto coarsed = coarsify(mp, 0); + if (isDead(mp, coarsed, mp_valid_per_bx)) + continue; + auto index = killTps(mp, coarsed, BX, mp_valid_per_bx); + mp_valid_per_bx[BX][index] = valid_cor_tp_t({true, mp, coarsed[3], coarsed[4], coarsed[5]}); + } + + std::vector outTPs; + for (auto &elem : mp_valid_per_bx) { + for (auto &mp_valid : elem.second) { + if (mp_valid.valid) { + outTPs.push_back(mp_valid.mp); + } + } + } + + for (auto &mp : SL2mps) + outTPs.push_back(mp); + return outTPs; +} + +std::vector MPCorFilter::coarsify(cmsdt::metaPrimitive mp, int sl) { + float pos_ch_f = mp.x; + + // translating into tdc counts + int pos_ch = int(round(pos_ch_f)); + int slope = (int)(mp.tanPhi); + + std::vector t0_slv, t0_coarse, pos_slv, pos_coarse, slope_slv, slope_coarse; + vhdl_int_to_unsigned(mp.t0, t0_slv); + vhdl_int_to_signed(pos_ch, pos_slv); + vhdl_int_to_signed(slope, slope_slv); + + vhdl_resize_unsigned(t0_slv, WIDTH_FULL_TIME); + vhdl_resize_signed(pos_slv, WIDTH_FULL_POS); + vhdl_resize_signed(slope_slv, WIDTH_FULL_SLOPE); + + t0_coarse = vhdl_slice(t0_slv, FSEG_T0_BX_LSB + 4, FSEG_T0_DISCARD_LSB); + pos_coarse = vhdl_slice(pos_slv, WIDTH_FULL_POS - 1, FSEG_POS_DISCARD_LSB); + slope_coarse = vhdl_slice(slope_slv, WIDTH_FULL_SLOPE - 1, FSEG_SLOPE_DISCARD_LSB); + + std::vector results; + int t0_coarse_int = vhdl_unsigned_to_int(t0_coarse); + int pos_coarse_int = vhdl_signed_to_int(pos_coarse); + int slope_coarse_int = vhdl_signed_to_int(slope_coarse); + + for (int index = 0; index <= 2; index++) { + auto aux_t0_coarse_int = + (t0_coarse_int + (index - 1)) % (int)std::pow(2, FSEG_T0_BX_LSB + 4 - (FSEG_T0_DISCARD_LSB)); + auto aux_pos_coarse_int = pos_coarse_int + (index - 1); + auto aux_slope_coarse_int = slope_coarse_int + (index - 1); + results.push_back(aux_t0_coarse_int); + results.push_back(aux_pos_coarse_int); + results.push_back(aux_slope_coarse_int); + } + return results; +} + +int MPCorFilter::match(cmsdt::metaPrimitive mp, std::vector coarsed, valid_cor_tp_t valid_cor_tp2) { + bool matched = ((coarsed[0] == valid_cor_tp2.coarsed_t0 || coarsed[3] == valid_cor_tp2.coarsed_t0 || + coarsed[6] == valid_cor_tp2.coarsed_t0) && + (coarsed[1] == valid_cor_tp2.coarsed_pos || coarsed[4] == valid_cor_tp2.coarsed_pos || + coarsed[7] == valid_cor_tp2.coarsed_pos) && + (coarsed[2] == valid_cor_tp2.coarsed_slope || coarsed[5] == valid_cor_tp2.coarsed_slope || + coarsed[8] == valid_cor_tp2.coarsed_slope)) && + (abs(mp.t0 / 25 - valid_cor_tp2.mp.t0 / 25) <= 1); + return ((int)matched) * 2 + (int)(mp.quality > valid_cor_tp2.mp.quality) + + (int)(mp.quality == valid_cor_tp2.mp.quality) * (int)(get_chi2(mp) < get_chi2(valid_cor_tp2.mp)); +} + +bool MPCorFilter::isDead(cmsdt::metaPrimitive mp, + std::vector coarsed, + std::map tps_per_bx) { + for (auto &elem : tps_per_bx) { + for (auto &mp_valid : elem.second) { + if (!mp_valid.valid) + continue; + int isMatched = match(mp, coarsed, mp_valid); + if (isMatched == 2) + return true; // matched and quality <= stored tp + } + } + return false; +} + +int MPCorFilter::killTps(cmsdt::metaPrimitive mp, + std::vector coarsed, + int bx, + std::map &tps_per_bx) { + int index_to_occupy = -1; + int index_to_kill = -1; + for (auto &elem : tps_per_bx) { + if (abs(bx - elem.first) > 2) + continue; + for (size_t i = 0; i < elem.second.size(); i++) { + if (elem.second[i].valid == 1) { + int isMatched = match(mp, coarsed, elem.second[i]); + if (isMatched == 3) { + elem.second[i].valid = false; + if (elem.first == bx && index_to_kill == -1) + index_to_kill = i; + } + } else if (elem.first == bx && index_to_occupy == -1) + index_to_occupy = i; + } + } + // My first option is to replace the one from my BX that I killed first + if (index_to_kill != -1) + return index_to_kill; + // If I wasn't able to kill anyone from my BX, I fill the first empty space + return index_to_occupy; +} + +int MPCorFilter::get_chi2(cmsdt::metaPrimitive mp) { + // chi2 is coarsified to the index of the chi2's highest bit set to 1 + + int chi2 = (int)round(mp.chi2 / (std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100)); + + std::vector chi2_unsigned, chi2_unsigned_msb; + vhdl_int_to_unsigned(chi2, chi2_unsigned); + + for (int i = (int)chi2_unsigned.size() - 1; i >= 0; i--) { + if (chi2_unsigned[i] == 1) { + return i; + } + } + return -1; +} + +void MPCorFilter::printmP(metaPrimitive mP) { + DTSuperLayerId slId(mP.rawId); + LogDebug("MPCorFilter") << slId << "\t" + << " " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 << " " << setw(2) + << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(5) << left << mP.tdc1 + << " " << setw(5) << left << mP.tdc2 << " " << setw(5) << left << mP.tdc3 << " " << setw(5) + << left << mP.tdc4 << " " << setw(10) << right << mP.x << " " << setw(9) << left << mP.tanPhi + << " " << setw(5) << left << mP.t0 << " " << setw(13) << left << mP.chi2; +} diff --git a/L1Trigger/DTTriggerPhase2/src/MPFilter.cc b/L1Trigger/DTTriggerPhase2/src/MPFilter.cc index c098972be4e89..fb45a0803f21c 100644 --- a/L1Trigger/DTTriggerPhase2/src/MPFilter.cc +++ b/L1Trigger/DTTriggerPhase2/src/MPFilter.cc @@ -8,6 +8,17 @@ using namespace std; // ============================================================================ MPFilter::MPFilter(const ParameterSet& pset) : debug_(pset.getUntrackedParameter("debug")) { // Obtention of parameters + int wh, st, se, maxdrift; + maxdrift_filename_ = pset.getParameter("maxdrift_filename"); + std::ifstream ifind(maxdrift_filename_.fullPath()); + if (ifind.fail()) { + throw cms::Exception("Missing Input File") + << "MPSLFilter::MPSLFilter() - Cannot find " << maxdrift_filename_.fullPath(); + } + while (ifind.good()) { + ifind >> wh >> st >> se >> maxdrift; + maxdriftinfo_[wh][st][se] = maxdrift; + } } MPFilter::~MPFilter() {} diff --git a/L1Trigger/DTTriggerPhase2/src/MPQualityEnhancerFilter.cc b/L1Trigger/DTTriggerPhase2/src/MPQualityEnhancerFilter.cc index 94b72556b9dc3..0a76fa7cffa6f 100644 --- a/L1Trigger/DTTriggerPhase2/src/MPQualityEnhancerFilter.cc +++ b/L1Trigger/DTTriggerPhase2/src/MPQualityEnhancerFilter.cc @@ -120,11 +120,7 @@ void MPQualityEnhancerFilter::filterCousins(std::vector &inMPaths bestI = i; } } - bool add_paths = (i == (int)(inMPaths.size() - 1)); - if (!add_paths) { - add_paths = areCousins(inMPaths[i], inMPaths[i + 1]) == 0; - } - if (!add_paths) { + if (areCousins(inMPaths[i], inMPaths[i + 1]) != 0) { primo_index++; } else { //areCousing==0 if (oneof4) { diff --git a/L1Trigger/DTTriggerPhase2/src/MPSLFilter.cc b/L1Trigger/DTTriggerPhase2/src/MPSLFilter.cc new file mode 100644 index 0000000000000..808eb38126d3f --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MPSLFilter.cc @@ -0,0 +1,278 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MPSLFilter.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +using namespace edm; +using namespace std; +using namespace cmsdt; + +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MPSLFilter::MPSLFilter(const ParameterSet &pset) : MPFilter(pset), debug_(pset.getUntrackedParameter("debug")) {} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void MPSLFilter::initialise(const edm::EventSetup &iEventSetup) {} + +void MPSLFilter::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector &inMPaths, + std::vector &outMPaths) { + if (debug_) + LogDebug("MPSLFilter") << "MPSLFilter: run"; + if (!inMPaths.empty()) { + int dum_sl_rawid = inMPaths[0].rawId; + DTSuperLayerId dumSlId(dum_sl_rawid); + DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector()); + max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1]; + DTSuperLayerId sl1Id(ChId.rawId(), 1); + DTSuperLayerId sl2Id(ChId.rawId(), 2); + DTSuperLayerId sl3Id(ChId.rawId(), 3); + + std::vector SL1metaPrimitives; + std::vector SL2metaPrimitives; + std::vector SL3metaPrimitives; + for (const auto &metaprimitiveIt : inMPaths) { + // int BX = metaprimitiveIt.t0 / 25; + if (metaprimitiveIt.rawId == sl1Id.rawId()) + SL1metaPrimitives.push_back(metaprimitiveIt); + else if (metaprimitiveIt.rawId == sl3Id.rawId()) + SL3metaPrimitives.push_back(metaprimitiveIt); + else if (metaprimitiveIt.rawId == sl2Id.rawId()) + SL2metaPrimitives.push_back(metaprimitiveIt); + } + + auto filteredSL1MPs = filter(SL1metaPrimitives); + auto filteredSL2MPs = filter(SL2metaPrimitives); + auto filteredSL3MPs = filter(SL3metaPrimitives); + + for (auto &mp : filteredSL1MPs) + outMPaths.push_back(mp); + for (auto &mp : filteredSL2MPs) + outMPaths.push_back(mp); + for (auto &mp : filteredSL3MPs) + outMPaths.push_back(mp); + } +} + +void MPSLFilter::finish(){}; + +/////////////////////////// +/// OTHER METHODS + +std::vector MPSLFilter::filter(std::vector mps) { + std::map mp_valid_per_bx; + for (auto &mp : mps) { + int BX = mp.t0 / 25; + if (mp_valid_per_bx.find(BX) == mp_valid_per_bx.end()) + mp_valid_per_bx[BX] = valid_tp_arr_t(6); + + // is this mp getting killed? + if (isDead(mp, mp_valid_per_bx)) + continue; + // if not, let's kill other mps + auto index = killTps(mp, BX, mp_valid_per_bx); + if (index == -1) + continue; + mp_valid_per_bx[BX][index] = valid_tp_t({true, mp}); + } + + std::vector outTPs; + for (auto &elem : mp_valid_per_bx) { + for (auto &mp_valid : elem.second) { + if (mp_valid.valid) + outTPs.push_back(mp_valid.mp); + } + } + + return outTPs; +} + +int MPSLFilter::match(cmsdt::metaPrimitive mp, cmsdt::metaPrimitive mp2) { + if ((mp.quality == mp2.quality) && (mp.quality == LOWQ || mp2.quality == CLOWQ)) + return 1; + + // CONFIRMATION, FIXME /////////////////////////// + // if (mp.quality == CLOWQ && mp2.quality == HIGHQ) { + // if (share_hit(mp, mp2)) return 2; + // return 3; + // } + // if (mp.quality == HIGHQ && mp2.quality == CLOWQ) { + // if (share_hit(mp, mp2)) return 4; + // return 5; + // } + ////////////////////////////////////////////////// + + if (mp.quality > mp2.quality) { + if (share_hit(mp, mp2)) + return 2; + return 3; + } + if (mp.quality < mp2.quality) { + if (share_hit(mp, mp2)) + return 4; + return 5; + } + if (share_hit(mp, mp2)) { + if (smaller_chi2(mp, mp2) == 0) + return 6; + return 7; + } + if (smaller_chi2(mp, mp2) == 0) + return 8; + return 9; +} + +bool MPSLFilter::isDead(cmsdt::metaPrimitive mp, std::map tps_per_bx) { + for (auto &elem : tps_per_bx) { + for (auto &mp_valid : elem.second) { + if (!mp_valid.valid) + continue; + int isMatched = match(mp, mp_valid.mp); + if (isMatched == 4 || isMatched == 7) + return true; + } + } + return false; +} + +int MPSLFilter::smaller_chi2(cmsdt::metaPrimitive mp, cmsdt::metaPrimitive mp2) { + auto chi2_1 = get_chi2(mp); + auto chi2_2 = get_chi2(mp2); + if (chi2_1 < chi2_2) + return 0; + return 1; +} + +int MPSLFilter::get_chi2(cmsdt::metaPrimitive mp) { + // CHI2 is converted to an unsigned in which 4 msb are the exponent + // of a float-like value and the rest of the bits are the mantissa + // (without the first 1). So comparing these reduced-width unsigned + // values is equivalent to comparing rounded versions of the chi2 + + int chi2 = (int)round(mp.chi2 / (std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100)); + + std::vector chi2_unsigned, chi2_unsigned_msb; + vhdl_int_to_unsigned(chi2, chi2_unsigned); + + if (chi2_unsigned.size() > 2) { + for (int i = (int)chi2_unsigned.size() - 1; i >= 2; i--) { + if (chi2_unsigned[i] == 1) { + vhdl_int_to_unsigned(i - 1, chi2_unsigned_msb); + + for (int j = i - 1; j > i - 3; j--) { + chi2_unsigned_msb.insert(chi2_unsigned_msb.begin(), chi2_unsigned[j]); + } + return vhdl_unsigned_to_int(chi2_unsigned_msb); + } + } + } + vhdl_resize_unsigned(chi2_unsigned, 2); + return vhdl_unsigned_to_int(vhdl_slice(chi2_unsigned, 1, 0)); +} + +int MPSLFilter::killTps(cmsdt::metaPrimitive mp, int bx, std::map &tps_per_bx) { + int index_to_occupy = -1; + int index_to_kill = -1; + for (auto &elem : tps_per_bx) { + if (abs(bx - elem.first) > 16) + continue; + for (size_t i = 0; i < elem.second.size(); i++) { + if (elem.second[i].valid == 1) { + int isMatched = match(mp, elem.second[i].mp); + if (isMatched == 2 || isMatched == 6) { + elem.second[i].valid = false; + if (elem.first == bx && index_to_kill == -1) + index_to_kill = i; + } + } else if (elem.first == bx && index_to_occupy == -1) + index_to_occupy = i; + } + } + // My first option is to replace the one from my BX that I killed first + if (index_to_kill != -1) + return index_to_kill; + // If I wasn't able to kill anyone from my BX, I fill the first empty space + if (index_to_occupy != -1) + return index_to_occupy; + // If I'm a 3h and there were no empty spaces, I don't replace any tp + if (mp.quality == LOWQ) + return -1; + // If I'm a 4h, I replace the first 3h or the 4h with the biggest chi2. + // Let's try to find both + int biggest_chi2 = 0; + int clowq_index = -1; + for (size_t i = 0; i < tps_per_bx[bx].size(); i++) { + if (tps_per_bx[bx][i].mp.quality == LOWQ) + return i; + if (tps_per_bx[bx][i].mp.quality == CLOWQ && clowq_index == -1) { + clowq_index = i; + continue; + } + auto chi2 = get_chi2(tps_per_bx[bx][i].mp); + if (chi2 > biggest_chi2) { + index_to_kill = i; + biggest_chi2 = chi2; + } + } + // If I found a confirmed 3h, I replace that one + if (clowq_index != -1) + return clowq_index; + // If all stored tps are 4h and their chi2 is smaller than mine, I don't replace any + if (biggest_chi2 < get_chi2(mp)) + return -1; + // If at least one chi2 is bigger than mine, I replace the corresponding tp + return index_to_kill; +} + +int MPSLFilter::share_hit(cmsdt::metaPrimitive mp, cmsdt::metaPrimitive mp2) { + // This function returns the layer % 4 (1 to 4) of the hit that is shared between TPs + // If they don't share any hits or the last hit of the latest one differs in more than + // SLFILT_MAX_SEG1T0_TO_SEG2ARRIVAL w.r.t. the t0 of the other, returns 0 + + // checking that they are from the same SL + if (mp.rawId != mp2.rawId) + return 0; + + bool isSL1 = ((int)(mp2.wi1 != -1) + (int)(mp2.wi2 != -1) + (int)(mp2.wi3 != -1) + (int)(mp2.wi4 != -1)) >= 3; + + int tdc_mp[NUM_LAYERS_2SL] = {mp.tdc1, mp.tdc2, mp.tdc3, mp.tdc4, mp.tdc5, mp.tdc6, mp.tdc7, mp.tdc8}; + int tdc_mp2[NUM_LAYERS_2SL] = {mp2.tdc1, mp2.tdc2, mp2.tdc3, mp2.tdc4, mp2.tdc5, mp2.tdc6, mp2.tdc7, mp2.tdc8}; + int max_tdc_mp = -999, max_tdc_mp2 = -999; + + for (size_t i = 0; i < NUM_LAYERS_2SL; i++) { + if (tdc_mp[i] > max_tdc_mp) + max_tdc_mp = tdc_mp[i]; + if (tdc_mp2[i] > max_tdc_mp2) + max_tdc_mp2 = tdc_mp2[i]; + } + + if (mp.t0 / LHC_CLK_FREQ + SLFILT_MAX_SEG1T0_TO_SEG2ARRIVAL < max_tdc_mp2 / LHC_CLK_FREQ || + mp2.t0 / LHC_CLK_FREQ + SLFILT_MAX_SEG1T0_TO_SEG2ARRIVAL < max_tdc_mp / LHC_CLK_FREQ) + return 0; + + if ((isSL1 && (mp.wi1 == mp2.wi1 and mp.tdc1 == mp2.tdc1 and mp.wi1 != -1 and mp.tdc1 != -1)) || + (!isSL1 && (mp.wi5 == mp2.wi5 and mp.tdc5 == mp2.tdc5 and mp.wi5 != -1 and mp.tdc5 != -1))) + return 1; + if ((isSL1 && (mp.wi2 == mp2.wi2 and mp.tdc2 == mp2.tdc2 and mp.wi2 != -1 and mp.tdc2 != -1)) || + (!isSL1 && (mp.wi6 == mp2.wi6 and mp.tdc6 == mp2.tdc6 and mp.wi6 != -1 and mp.tdc6 != -1))) + return 2; + if ((isSL1 && (mp.wi3 == mp2.wi3 and mp.tdc3 == mp2.tdc3 and mp.wi3 != -1 and mp.tdc3 != -1)) || + (!isSL1 && (mp.wi7 == mp2.wi7 and mp.tdc7 == mp2.tdc7 and mp.wi7 != -1 and mp.tdc7 != -1))) + return 3; + if ((isSL1 && (mp.wi4 == mp2.wi4 and mp.tdc4 == mp2.tdc4 and mp.wi4 != -1 and mp.tdc4 != -1)) || + (!isSL1 && (mp.wi8 == mp2.wi8 and mp.tdc8 == mp2.tdc8 and mp.wi8 != -1 and mp.tdc8 != -1))) + return 4; + return 0; +} + +void MPSLFilter::printmP(metaPrimitive mP) { + DTSuperLayerId slId(mP.rawId); + LogDebug("MPSLFilter") << slId << "\t" + << " " << setw(2) << left << mP.wi1 << " " << setw(2) << left << mP.wi2 << " " << setw(2) + << left << mP.wi3 << " " << setw(2) << left << mP.wi4 << " " << setw(5) << left << mP.tdc1 + << " " << setw(5) << left << mP.tdc2 << " " << setw(5) << left << mP.tdc3 << " " << setw(5) + << left << mP.tdc4 << " " << setw(10) << right << mP.x << " " << setw(9) << left << mP.tanPhi + << " " << setw(5) << left << mP.t0 << " " << setw(13) << left << mP.chi2; +} diff --git a/L1Trigger/DTTriggerPhase2/src/MuonPathAssociator.cc b/L1Trigger/DTTriggerPhase2/src/MuonPathAssociator.cc index e96f1634e2980..0f406b32c47bf 100644 --- a/L1Trigger/DTTriggerPhase2/src/MuonPathAssociator.cc +++ b/L1Trigger/DTTriggerPhase2/src/MuonPathAssociator.cc @@ -108,7 +108,7 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, SL3metaPrimitives.push_back(metaprimitiveIt); } - if (SL1metaPrimitives.empty() or SL3metaPrimitives.empty()) + if (SL1metaPrimitives.empty() and SL3metaPrimitives.empty()) continue; if (debug_) @@ -119,12 +119,12 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, bool at_least_one_SL1_confirmation = false; bool at_least_one_SL3_confirmation = false; - vector useFitSL1; + bool useFitSL1[SL1metaPrimitives.size()]; for (unsigned int i = 0; i < SL1metaPrimitives.size(); i++) - useFitSL1.push_back(false); - vector useFitSL3; + useFitSL1[i] = false; + bool useFitSL3[SL3metaPrimitives.size()]; for (unsigned int i = 0; i < SL3metaPrimitives.size(); i++) - useFitSL3.push_back(false); + useFitSL3[i] = false; //SL1-SL3 vector chamberMetaPrimitives; @@ -133,7 +133,7 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, int sl1 = 0; int sl3 = 0; for (auto SL1metaPrimitive = SL1metaPrimitives.begin(); SL1metaPrimitive != SL1metaPrimitives.end(); - ++SL1metaPrimitive, sl1++, sl3 = 0) { + ++SL1metaPrimitive, sl1++, sl3 = -1) { if (clean_chi2_correlation_) at_least_one_correlation = false; for (auto SL3metaPrimitive = SL3metaPrimitives.begin(); SL3metaPrimitive != SL3metaPrimitives.end(); @@ -412,13 +412,19 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, next_tdc = best_tdc; next_layer = best_layer; next_lat = best_lat; + + best_wire = (*digiIt).wire(); + best_tdc = (*digiIt).time(); + best_layer = dtLId.layer(); + best_lat = lat; matched_digis++; + } else if (dtLId.layer() == + best_layer) { // same layer than stored, just substituting the hit, no matched_digis++; + best_wire = (*digiIt).wire(); + best_tdc = (*digiIt).time(); + best_layer = dtLId.layer(); + best_lat = lat; } - best_wire = (*digiIt).wire(); - best_tdc = (*digiIt).time(); - best_layer = dtLId.layer(); - best_lat = lat; - } else if ((std::abs(x_inSL3 - x_wire) >= minx) && (std::abs(x_inSL3 - x_wire) < min2x)) { // same layer than the stored in best, no hit added if (dtLId.layer() == best_layer) @@ -427,7 +433,8 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, // buggy, as we could have stored as next LayerX -> LayerY -> LayerX, and this should // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more // makes no difference - matched_digis++; + else if (dtLId.layer() != next_layer) + matched_digis++; // whether the layer is the same for this hit and the stored in next, we substitute // the one stored and modify the min distance min2x = std::abs(x_inSL3 - x_wire); @@ -640,12 +647,19 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, next_tdc = best_tdc; next_layer = best_layer; next_lat = best_lat; + + best_wire = (*digiIt).wire(); + best_tdc = (*digiIt).time(); + best_layer = dtLId.layer(); + best_lat = lat; matched_digis++; + } else if (dtLId.layer() == + best_layer) { // same layer than stored, just substituting the hit, no matched_digis++; + best_wire = (*digiIt).wire(); + best_tdc = (*digiIt).time(); + best_layer = dtLId.layer(); + best_lat = lat; } - best_wire = (*digiIt).wire(); - best_tdc = (*digiIt).time(); - best_layer = dtLId.layer(); - best_lat = lat; } else if ((std::abs(x_inSL1 - x_wire) >= minx) && (std::abs(x_inSL1 - x_wire) < min2x)) { // same layer than the stored in best, no hit added if (dtLId.layer() == best_layer) @@ -654,7 +668,8 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, // buggy, as we could have stored as next LayerX -> LayerY -> LayerX, and this should // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more // makes no difference - matched_digis++; + else if (dtLId.layer() != next_layer) + matched_digis++; // whether the layer is the same for this hit and the stored in next, we substitute // the one stored and modify the min distance min2x = std::abs(x_inSL1 - x_wire); @@ -986,9 +1001,9 @@ void MuonPathAssociator::correlateMPaths(edm::Handle dtdigis, } void MuonPathAssociator::removeSharingFits(vector &chamberMPaths, vector &allMPaths) { - vector useFit; + bool useFit[chamberMPaths.size()]; for (unsigned int i = 0; i < chamberMPaths.size(); i++) { - useFit.push_back(true); + useFit[i] = true; } for (unsigned int i = 0; i < chamberMPaths.size(); i++) { if (debug_) diff --git a/L1Trigger/DTTriggerPhase2/src/MuonPathConfirmator.cc b/L1Trigger/DTTriggerPhase2/src/MuonPathConfirmator.cc new file mode 100644 index 0000000000000..8390b62fe8632 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MuonPathConfirmator.cc @@ -0,0 +1,237 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; + +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MuonPathConfirmator::MuonPathConfirmator(const ParameterSet &pset, edm::ConsumesCollector &iC) + : debug_(pset.getUntrackedParameter("debug")), + minx_match_2digis_(pset.getParameter("minx_match_2digis")) { + if (debug_) + LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: constructor"; + + //shift phi + int rawId; + shift_filename_ = pset.getParameter("shift_filename"); + std::ifstream ifin3(shift_filename_.fullPath()); + double shift; + if (ifin3.fail()) { + throw cms::Exception("Missing Input File") + << "MuonPathConfirmator::MuonPathConfirmator() - Cannot find " << shift_filename_.fullPath(); + } + while (ifin3.good()) { + ifin3 >> rawId >> shift; + shiftinfo_[rawId] = shift; + } + + int wh, st, se, maxdrift; + maxdrift_filename_ = pset.getParameter("maxdrift_filename"); + std::ifstream ifind(maxdrift_filename_.fullPath()); + if (ifind.fail()) { + throw cms::Exception("Missing Input File") + << "MPSLFilter::MPSLFilter() - Cannot find " << maxdrift_filename_.fullPath(); + } + while (ifind.good()) { + ifind >> wh >> st >> se >> maxdrift; + maxdriftinfo_[wh][st][se] = maxdrift; + } +} + +MuonPathConfirmator::~MuonPathConfirmator() { + if (debug_) + LogDebug("MuonPathConfirmator") << "MuonPathAnalyzer: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ + +void MuonPathConfirmator::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + std::vector inMetaPrimitives, + edm::Handle dtdigis, + std::vector &outMetaPrimitives) { + if (debug_) + LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: run"; + + // fit per SL (need to allow for multiple outputs for a single mpath) + if (!inMetaPrimitives.empty()) { + int dum_sl_rawid = inMetaPrimitives[0].rawId; + DTSuperLayerId dumSlId(dum_sl_rawid); + DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector()); + max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1]; + } + + for (auto &mp : inMetaPrimitives) { + analyze(mp, dtdigis, outMetaPrimitives); + } +} + +void MuonPathConfirmator::initialise(const edm::EventSetup &iEventSetup) { + if (debug_) + LogDebug("MuonPathConfirmator") << "MuonPathConfirmator::initialiase"; +} + +void MuonPathConfirmator::finish() { + if (debug_) + LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: finish"; +}; + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +void MuonPathConfirmator::analyze(cmsdt::metaPrimitive mp, + edm::Handle dtdigis, + std::vector &outMetaPrimitives) { + int dum_sl_rawid = mp.rawId; + DTSuperLayerId dumSlId(dum_sl_rawid); + DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector()); + DTSuperLayerId sl1Id(ChId.rawId(), 1); + DTSuperLayerId sl3Id(ChId.rawId(), 3); + + DTWireId wireIdSL1(sl1Id, 2, 1); + DTWireId wireIdSL3(sl3Id, 2, 1); + auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()]; + bool isSL1 = (mp.rawId == sl1Id.rawId()); + bool isSL3 = (mp.rawId == sl3Id.rawId()); + if (!isSL1 && !isSL3) + outMetaPrimitives.emplace_back(mp); + else { + int best_tdc = -1; + int next_tdc = -1; + int best_wire = -1; + int next_wire = -1; + int best_layer = -1; + int next_layer = -1; + int best_lat = -1; + int next_lat = -1; + int lat = -1; + int matched_digis = 0; + + int position_prec = ((int)(mp.x)) << PARTIALS_PRECISSION; + int slope_prec = ((int)(mp.tanPhi)) << PARTIALS_PRECISSION; + + int slope_x_halfchamb = (((long int)slope_prec) * SEMICHAMBER_H) >> SEMICHAMBER_RES_SHR; + int slope_x_3semicells = (slope_prec * 3) >> LYRANDAHALF_RES_SHR; + int slope_x_1semicell = (slope_prec * 1) >> LYRANDAHALF_RES_SHR; + + for (const auto &dtLayerId_It : *dtdigis) { + const DTLayerId dtLId = dtLayerId_It.first; + // creating a new DTSuperLayerId object to compare with the required SL id + const DTSuperLayerId dtSLId(dtLId.wheel(), dtLId.station(), dtLId.sector(), dtLId.superLayer()); + bool hitFromSL1 = (dtSLId.rawId() == sl1Id.rawId()); + bool hitFromSL3 = (dtSLId.rawId() == sl3Id.rawId()); + if (!(hitFromSL1 || hitFromSL3)) // checking hits are from one of the other SL of the same chamber + continue; + double minx = 10 * minx_match_2digis_ * ((double)max_drift_tdc / (double)CELL_SEMILENGTH); + double min2x = 10 * minx_match_2digis_ * ((double)max_drift_tdc / (double)CELL_SEMILENGTH); + if (isSL1 != hitFromSL1) { // checking hits have the opposite SL than the TP + for (auto digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; ++digiIt) { + if ((*digiIt).time() < mp.t0) + continue; + int wp_semicells = ((*digiIt).wire() - 1 - SL1_CELLS_OFFSET) * 2 + 1; + int ly = dtLId.layer() - 1; + if (ly % 2 == 1) + wp_semicells -= 1; + if (hitFromSL3) + wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH); + double hit_position = wp_semicells * max_drift_tdc + + ((*digiIt).time() - mp.t0) * (double)TIME_TO_TDC_COUNTS / (double)LHC_CLK_FREQ; + double hit_position_left = wp_semicells * max_drift_tdc - + ((*digiIt).time() - mp.t0) * (double)TIME_TO_TDC_COUNTS / (double)LHC_CLK_FREQ; + // extrapolating position to the layer of the hit + // mp.position is referred to the center between SLs, so one has to add half the distance between SLs + // + half a cell height to get to the first wire + ly * cell height to reach the desired ly + // 10 * VERT_PHI1_PHI3 / 2 + (CELL_HEIGHT / 2) + ly * CELL_HEIGHT = (10 * VERT_PHI1_PHI3 + (2 * ly + 1) * CELL_HEIGHT) / 2 + + int position_in_layer = position_prec + (1 - 2 * (int)hitFromSL1) * slope_x_halfchamb; + if (ly == 0) + position_in_layer -= slope_x_3semicells; + if (ly == 1) + position_in_layer -= slope_x_1semicell; + if (ly == 2) + position_in_layer += slope_x_1semicell; + if (ly == 3) + position_in_layer += slope_x_3semicells; + position_in_layer = position_in_layer >> PARTIALS_PRECISSION; + + if (std::abs(position_in_layer - hit_position_left) < std::abs(position_in_layer - hit_position)) { + lat = 0; + hit_position = hit_position_left; + } + if (std::abs(position_in_layer - hit_position) < minx) { + // different layer than the stored in best, hit added, matched_digis++;. This approach in somewhat + // buggy, as we could have stored as best LayerX -> LayerY -> LayerX, and this should + // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more + // makes no difference + if (dtLId.layer() != best_layer) { + minx = std::abs(position_in_layer - hit_position); + next_wire = best_wire; + next_tdc = best_tdc; + next_layer = best_layer; + next_lat = best_lat; + matched_digis++; + } + best_wire = (*digiIt).wire() - 1; + best_tdc = (*digiIt).time(); + best_layer = dtLId.layer(); + best_lat = lat; + + } else if ((std::abs(position_in_layer - hit_position) >= minx) && + (std::abs(position_in_layer - hit_position) < min2x)) { + // same layer than the stored in best, no hit added + if (dtLId.layer() == best_layer) + continue; + // different layer than the stored in next, hit added. This approach in somewhat + // buggy, as we could have stored as next LayerX -> LayerY -> LayerX, and this should + // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more + // makes no difference + matched_digis++; + // whether the layer is the same for this hit and the stored in next, we substitute + // the one stored and modify the min distance + min2x = std::abs(position_in_layer - hit_position); + next_wire = (*digiIt).wire() - 1; + next_tdc = (*digiIt).time(); + next_layer = dtLId.layer(); + next_lat = lat; + } + } + } + } + int new_quality = mp.quality; + std::vector wi_c(4, -1), tdc_c(4, -1), lat_c(4, -1); + if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) { // actually confirm + new_quality = CHIGHQ; + if (mp.quality == LOWQ) + new_quality = CLOWQ; + + wi_c[next_layer - 1] = next_wire; + tdc_c[next_layer - 1] = next_tdc; + lat_c[next_layer - 1] = next_lat; + + wi_c[best_layer - 1] = best_wire; + tdc_c[best_layer - 1] = best_tdc; + lat_c[best_layer - 1] = best_lat; + } + if (isSL1) { + outMetaPrimitives.emplace_back(metaPrimitive( + {mp.rawId, mp.t0, mp.x, mp.tanPhi, mp.phi, mp.phiB, mp.phi_cmssw, mp.phiB_cmssw, mp.chi2, new_quality, + mp.wi1, mp.tdc1, mp.lat1, mp.wi2, mp.tdc2, mp.lat2, mp.wi3, mp.tdc3, mp.lat3, mp.wi4, + mp.tdc4, mp.lat4, wi_c[0], tdc_c[0], lat_c[0], wi_c[1], tdc_c[1], lat_c[1], wi_c[2], tdc_c[2], + lat_c[2], wi_c[3], tdc_c[3], lat_c[3], -1})); + } else { + outMetaPrimitives.emplace_back( + metaPrimitive({mp.rawId, mp.t0, mp.x, mp.tanPhi, mp.phi, mp.phiB, mp.phi_cmssw, + mp.phiB_cmssw, mp.chi2, new_quality, wi_c[0], tdc_c[0], lat_c[0], wi_c[1], + tdc_c[1], lat_c[1], wi_c[2], tdc_c[2], lat_c[2], wi_c[3], tdc_c[3], + lat_c[3], mp.wi5, mp.tdc5, mp.lat5, mp.wi6, mp.tdc6, mp.lat6, + mp.wi7, mp.tdc7, mp.lat7, mp.wi8, mp.tdc8, mp.lat8, -1})); + } + } //SL2 +} diff --git a/L1Trigger/DTTriggerPhase2/src/MuonPathCorFitter.cc b/L1Trigger/DTTriggerPhase2/src/MuonPathCorFitter.cc new file mode 100644 index 0000000000000..4ef798f843119 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MuonPathCorFitter.cc @@ -0,0 +1,490 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathCorFitter.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MuonPathCorFitter::MuonPathCorFitter(const ParameterSet& pset, + edm::ConsumesCollector& iC, + std::shared_ptr& globalcoordsobtainer) + : MuonPathFitter(pset, iC, globalcoordsobtainer), dT0_correlate_TP_(pset.getParameter("dT0_correlate_TP")) { + if (debug_) + LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: constructor"; + + // LUTs + both_sl_filename_ = pset.getParameter("lut_2sl"); + + fillLuts(); + + setChi2Th(pset.getParameter("chi2corTh")); + setTanPhiTh(pset.getParameter("dTanPsi_correlate_TP")); +} + +MuonPathCorFitter::~MuonPathCorFitter() { + if (debug_) + LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void MuonPathCorFitter::initialise(const edm::EventSetup& iEventSetup) { + if (debug_) + LogDebug("MuonPathCorFitter") << "MuonPathCorFitter::initialiase"; + + auto geom = iEventSetup.getHandle(dtGeomH); + dtGeo_ = &(*geom); +} + +void MuonPathCorFitter::run(edm::Event& iEvent, + const edm::EventSetup& iEventSetup, + std::vector& inMPaths, + std::vector& outMPaths) { + if (debug_) + LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: run"; + if (!inMPaths.empty()) { + int dum_sl_rawid = inMPaths[0].rawId; + DTSuperLayerId dumSlId(dum_sl_rawid); + DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector()); + max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1]; + DTSuperLayerId sl1Id(ChId.rawId(), 1); + DTSuperLayerId sl3Id(ChId.rawId(), 3); + + std::map> SL1metaPrimitivesPerBX; + std::map> SL3metaPrimitivesPerBX; + for (const auto& metaprimitiveIt : inMPaths) { + int BX = metaprimitiveIt.t0 / 25; + if (metaprimitiveIt.rawId == sl1Id.rawId()) + SL1metaPrimitivesPerBX[BX].push_back(metaprimitiveIt); + else if (metaprimitiveIt.rawId == sl3Id.rawId()) + SL3metaPrimitivesPerBX[BX].push_back(metaprimitiveIt); + } + + std::vector bxs_to_consider; + bxs_to_consider.reserve(SL1metaPrimitivesPerBX.size()); + for (auto& prims_sl1 : SL1metaPrimitivesPerBX) + bxs_to_consider.push_back(bx_sl_vector({prims_sl1.first, prims_sl1.second, 1})); + + for (auto& prims_sl3 : SL3metaPrimitivesPerBX) + bxs_to_consider.push_back(bx_sl_vector({prims_sl3.first, prims_sl3.second, 3})); + + std::stable_sort(bxs_to_consider.begin(), bxs_to_consider.end(), bxSort); + + std::vector mps_q8; + std::vector mps_q7; + std::vector mps_q6; + + for (size_t ibx = 1; ibx < bxs_to_consider.size(); ibx++) { + for (size_t ibx2 = 0; ibx2 < ibx; ibx2++) { + if (bxs_to_consider[ibx].sl != bxs_to_consider[ibx2].sl && + (abs(bxs_to_consider[ibx].bx - bxs_to_consider[ibx2].bx)) <= MAX_BX_FOR_COR) { + int isl1 = 0; + for (auto& prim1 : bxs_to_consider[ibx].mps) { + if (isl1 >= MAX_PRIM_PER_BX_FOR_COR) + break; + int isl2 = 0; + for (auto& prim2 : bxs_to_consider[ibx2].mps) { + if (isl2 >= MAX_PRIM_PER_BX_FOR_COR) + break; + if (bxs_to_consider[ibx].sl == 1) { + if (!canCorrelate(prim1, prim2)) { + continue; + } + if (prim1.quality >= 3 && prim2.quality >= 3) + mps_q8.push_back(mp_group({prim1, prim2})); + else if ((prim1.quality >= 3 && prim2.quality < 3) || (prim1.quality < 3 && prim2.quality >= 3)) + mps_q7.push_back(mp_group({prim1, prim2})); + else + mps_q6.push_back(mp_group({prim1, prim2})); + } else { + if (!canCorrelate(prim2, prim1)) { + continue; + } + if (prim2.quality >= 3 && prim1.quality >= 3) + mps_q8.push_back(mp_group({prim2, prim1})); + else if ((prim2.quality >= 3 && prim1.quality < 3) || (prim2.quality < 3 && prim1.quality >= 3)) + mps_q7.push_back(mp_group({prim2, prim1})); + else + mps_q6.push_back(mp_group({prim2, prim1})); + } + isl2++; + } + isl1++; + } + } + } // looping over the 0 -> N-1 BX groups + } // looping over the 1 -> N BX groups + int iq = 0; + for (size_t i = 0; i < mps_q8.size(); i++) { + if (iq >= MAX_PRIM_FOR_COR) + break; + analyze(mps_q8[i], outMPaths); + iq += 1; + } + for (size_t i = 0; i < mps_q7.size(); i++) { + if (iq >= MAX_PRIM_FOR_COR) + break; + analyze(mps_q7[i], outMPaths); + iq += 1; + } + for (size_t i = 0; i < mps_q6.size(); i++) { + if (iq >= MAX_PRIM_FOR_COR) + break; + analyze(mps_q6[i], outMPaths); + iq += 1; + } + } +} + +bool MuonPathCorFitter::canCorrelate(cmsdt::metaPrimitive mp_sl1, cmsdt::metaPrimitive mp_sl3) { + // moving position from SL RF to chamber RF + float pos_ch_sl1_f = mp_sl1.x; + float pos_ch_sl3_f = mp_sl3.x; + + // translating into tdc counts + int pos_ch_sl1 = int(pos_ch_sl1_f); + int pos_ch_sl3 = int(pos_ch_sl3_f); + + int slope_sl1 = (int)mp_sl1.tanPhi; + int slope_sl3 = (int)mp_sl3.tanPhi; + + if (abs((slope_sl1 >> WIDTH_POS_SLOPE_CORR) - (slope_sl3 >> WIDTH_POS_SLOPE_CORR)) > 1) + return false; + + if (abs((pos_ch_sl1 >> WIDTH_POS_SLOPE_CORR) - (pos_ch_sl3 >> WIDTH_POS_SLOPE_CORR)) > 1) + return false; + + if (abs(mp_sl1.t0 - mp_sl3.t0) > dT0_correlate_TP_) + return false; + + return true; +} + +void MuonPathCorFitter::finish() { + if (debug_) + LogDebug("MuonPathCorFitter") << "MuonPathCorFitter: finish"; +}; + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +void MuonPathCorFitter::analyze(mp_group mp, std::vector& metaPrimitives) { + //FIXME + DTSuperLayerId MuonPathSLId(mp[0].rawId); // SL1 + + DTChamberId ChId(MuonPathSLId.wheel(), MuonPathSLId.station(), MuonPathSLId.sector()); + + DTSuperLayerId MuonPathSL1Id(ChId.wheel(), ChId.station(), ChId.sector(), 1); + DTSuperLayerId MuonPathSL3Id(ChId.wheel(), ChId.station(), ChId.sector(), 3); + DTWireId wireIdSL1(MuonPathSL1Id, 2, 1); + DTWireId wireIdSL3(MuonPathSL3Id, 2, 1); + auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()]; + + fit_common_in_t fit_common_in; + + // 8-element vectors, for the 8 layers. As here we are fitting one SL only, we leave the other SL values as dummy ones + fit_common_in.hits = {}; + fit_common_in.hits_valid = {}; + short quality = 0; + if (mp[0].quality >= 3 && mp[1].quality >= 3) + quality = 8; + else if ((mp[0].quality >= 3 && mp[1].quality < 3) || (mp[0].quality < 3 && mp[1].quality >= 3)) + quality = 7; + else + quality = 6; + + std::vector missing_layers; + + for (int isl = 0; isl < 2; isl++) { + int wire[4], tdc[4]; + if (isl != 1) { + wire[0] = mp[isl].wi1; + tdc[0] = mp[isl].tdc1; + wire[1] = mp[isl].wi2; + tdc[1] = mp[isl].tdc2; + wire[2] = mp[isl].wi3; + tdc[2] = mp[isl].tdc3; + wire[3] = mp[isl].wi4; + tdc[3] = mp[isl].tdc4; + } else { + wire[0] = mp[isl].wi5; + tdc[0] = mp[isl].tdc5; + wire[1] = mp[isl].wi6; + tdc[1] = mp[isl].tdc6; + wire[2] = mp[isl].wi7; + tdc[2] = mp[isl].tdc7; + wire[3] = mp[isl].wi8; + tdc[3] = mp[isl].tdc8; + } + + for (int i = 0; i < NUM_LAYERS; i++) { + if (wire[i] != -1) { + // Include both valid and non-valid hits. Non-valid values can be whatever, leaving all as -1 to make debugging easier. + auto ti = tdc[i]; + if (ti != -1) + ti = (int)round(((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ) * ti); + auto wi = wire[i]; + auto ly = i; + + int wp_semicells = (wi - SL1_CELLS_OFFSET) * 2 + 1; + if (ly % 2 == 1) + wp_semicells -= 1; + if (isl == 1) // SL3 + wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH); + float wp_tdc = wp_semicells * max_drift_tdc; + int wp = (int)((long int)(round(wp_tdc * std::pow(2, WIREPOS_WIDTH))) / (int)std::pow(2, WIREPOS_WIDTH)); + + // wp in tdc counts (still in floating point) + fit_common_in.hits.push_back({ti, wi, ly, wp}); + // fill valids as well + fit_common_in.hits_valid.push_back(1); + } else { + missing_layers.push_back(isl * NUM_LAYERS + i); + fit_common_in.hits.push_back({-1, -1, -1, -1}); + fit_common_in.hits_valid.push_back(0); + } + } + } + + int smallest_time = 999999, tmp_coarse_wirepos_1 = -1, tmp_coarse_wirepos_3 = -1; + // coarse_bctr is the 12 MSB of the smallest tdc + for (int isl = 0; isl < 2; isl++) { + for (size_t i = 0; i < NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[NUM_LAYERS * isl + i] == 0) + continue; + else if (fit_common_in.hits[NUM_LAYERS * isl + i].ti < smallest_time) + smallest_time = fit_common_in.hits[NUM_LAYERS * isl + i].ti; + } + } + if (fit_common_in.hits_valid[NUM_LAYERS * 0 + 0] == 1) + tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * 0 + 0].wp; + else + tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * 0 + 1].wp; + if (fit_common_in.hits_valid[NUM_LAYERS * 1 + 3] == 1) + tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * 1 + 3].wp; + else + tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * 1 + 2].wp; + + tmp_coarse_wirepos_1 = tmp_coarse_wirepos_1 >> WIREPOS_NORM_LSB_IGNORED; + tmp_coarse_wirepos_3 = tmp_coarse_wirepos_3 >> WIREPOS_NORM_LSB_IGNORED; + + fit_common_in.coarse_bctr = smallest_time >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME); + fit_common_in.coarse_wirepos = (tmp_coarse_wirepos_1 + tmp_coarse_wirepos_3) >> 1; + + fit_common_in.lateralities.clear(); + + auto rom_addr = get_rom_addr(mp, missing_layers); + + coeffs_t coeffs = + RomDataConvert(lut_2sl[rom_addr], COEFF_WIDTH_COR_T0, COEFF_WIDTH_COR_POSITION, COEFF_WIDTH_COR_SLOPE, 0, 7); + + // Filling lateralities + for (int isl = 0; isl < 2; isl++) { + int lat[4]; + if (isl != 1) { + lat[0] = mp[isl].lat1; + lat[1] = mp[isl].lat2; + lat[2] = mp[isl].lat3; + lat[3] = mp[isl].lat4; + } else { + lat[0] = mp[isl].lat5; + lat[1] = mp[isl].lat6; + lat[2] = mp[isl].lat7; + lat[3] = mp[isl].lat8; + } + + for (size_t i = 0; i < NUM_LAYERS; i++) { + fit_common_in.lateralities.push_back(lat[i]); + } + } + + fit_common_in.coeffs = coeffs; + + auto fit_common_out = fit(fit_common_in, + XI_COR_WIDTH, + COEFF_WIDTH_COR_T0, + COEFF_WIDTH_COR_POSITION, + COEFF_WIDTH_COR_SLOPE, + PRECISSION_COR_T0, + PRECISSION_COR_POSITION, + PRECISSION_COR_SLOPE, + PROD_RESIZE_COR_T0, + PROD_RESIZE_COR_POSITION, + PROD_RESIZE_COR_SLOPE, + max_drift_tdc, + 0); + + if (fit_common_out.valid_fit == 1) { + float t0_f = ((float)fit_common_out.t0) * (float)LHC_CLK_FREQ / (float)TIME_TO_TDC_COUNTS; + float slope_f = -fit_common_out.slope * ((float)CELL_SEMILENGTH / max_drift_tdc) * (1) / (CELL_SEMIHEIGHT * 16.); + if (std::abs(slope_f) > tanPhiTh_) + return; + + DTWireId wireId(MuonPathSLId, 2, 1); + float pos_ch_f = (float)(fit_common_out.position) * ((float)CELL_SEMILENGTH / (float)max_drift_tdc) / 10; + pos_ch_f += (SL1_CELLS_OFFSET * CELL_LENGTH) / 10.; + pos_ch_f += shiftinfo_[wireId.rawId()]; + + float chi2_f = fit_common_out.chi2 * std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100; + + // obtention of global coordinates using luts + int pos = (int)(10 * (pos_ch_f - shiftinfo_[wireId.rawId()]) * INCREASED_RES_POS_POW); + int slope = (int)(-slope_f * INCREASED_RES_SLOPE_POW); + auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), 0, pos, slope); + float phi = global_coords[0]; + float phiB = global_coords[1]; + + // obtention of global coordinates using cmssw geometry + double z = 0; + if (ChId.station() == 3 or ChId.station() == 4) { + z += Z_SHIFT_MB4; + } + GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(pos_ch_f, 0., z)); + int thisec = ChId.sector(); + if (thisec == 13) + thisec = 4; + if (thisec == 14) + thisec = 10; + float phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1); + float psi = atan(slope_f); + float phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw; + metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(), + t0_f, + (double)fit_common_out.position, + (double)fit_common_out.slope, + phi, + phiB, + phi_cmssw, + phiB_cmssw, + chi2_f, + quality, + mp[0].wi1, + mp[0].tdc1, + mp[0].lat1, + mp[0].wi2, + mp[0].tdc2, + mp[0].lat2, + mp[0].wi3, + mp[0].tdc3, + mp[0].lat3, + mp[0].wi4, + mp[0].tdc4, + mp[0].lat4, + mp[1].wi5, + mp[1].tdc5, + mp[1].lat5, + mp[1].wi6, + mp[1].tdc6, + mp[1].lat6, + mp[1].wi7, + mp[1].tdc7, + mp[1].lat7, + mp[1].wi8, + mp[1].tdc8, + mp[1].lat8, + -1})); + } + return; +} + +void MuonPathCorFitter::fillLuts() { + std::ifstream ifin2sl(both_sl_filename_.fullPath()); + std::string line; + while (ifin2sl.good()) { + ifin2sl >> line; + + std::vector myNumbers; + for (size_t i = 0; i < line.size(); i++) { + // This converts the char into an int and pushes it into vec + myNumbers.push_back(line[i] - '0'); // The digits will be in the same order as before + } + std::reverse(myNumbers.begin(), myNumbers.end()); + lut_2sl.push_back(myNumbers); + } + + return; +} + +int MuonPathCorFitter::get_rom_addr(mp_group mps, std::vector missing_hits) { + std::vector lats = { + mps[0].lat1, mps[0].lat2, mps[0].lat3, mps[0].lat4, mps[1].lat5, mps[1].lat6, mps[1].lat7, mps[1].lat8}; + + std::vector rom_addr; + if (missing_hits.size() == 1) + rom_addr.push_back(1); + else + rom_addr.push_back(0); + + if (missing_hits.size() == 1) { // 7 layers fit + if (missing_hits[0] < 4) + rom_addr.push_back(0); // First SL has 4 hits (1) or 3 (0) + else + rom_addr.push_back(1); + if (missing_hits[0] % 4 == 0) { + rom_addr.push_back(0); + rom_addr.push_back(0); + } else if (missing_hits[0] % 4 == 1) { + rom_addr.push_back(0); + rom_addr.push_back(1); + } else if (missing_hits[0] % 4 == 2) { + rom_addr.push_back(1); + rom_addr.push_back(0); + } else { // missing_hits[0] == 3 + rom_addr.push_back(1); + rom_addr.push_back(1); + } + for (size_t ilat = 0; ilat < lats.size(); ilat++) { + if ((int)ilat == missing_hits[0]) // only applies to 3-hit, as in 4-hit missL=-1 + continue; + auto lat = lats[ilat]; + if (lat == -1) + lat = 0; + rom_addr.push_back(lat); + } + + } else if (missing_hits.empty()) { // 8 layers fit + for (size_t ilat = 0; ilat < lats.size(); ilat++) { + auto lat = lats[ilat]; + if (lat == -1) + lat = 0; + rom_addr.push_back(lat); + } + auto lat = lats[NUM_LAYERS + 3]; + if (lat == -1) + lat = 0; + rom_addr.push_back(lat); + rom_addr.push_back(lat); + + } else { // 6 layers fit + for (int i = missing_hits.size() - 1; i >= 0; i--) { + if (missing_hits[i] % 4 == 0) { + rom_addr.push_back(0); + rom_addr.push_back(0); + } else if (missing_hits[i] % 4 == 1) { + rom_addr.push_back(0); + rom_addr.push_back(1); + } else if (missing_hits[i] % 4 == 2) { + rom_addr.push_back(1); + rom_addr.push_back(0); + } else { // missing_hits[i] % 4 == 3 + rom_addr.push_back(1); + rom_addr.push_back(1); + } + } + for (size_t ilat = 0; ilat < lats.size(); ilat++) { + if ((int)ilat == missing_hits[0] || (int)ilat == missing_hits[1]) // only applies to 3-hit, as in 4-hit missL=-1 + continue; + auto lat = lats[ilat]; + if (lat == -1) + lat = 0; + rom_addr.push_back(lat); + } + } + std::reverse(rom_addr.begin(), rom_addr.end()); + return vhdl_unsigned_to_int(rom_addr); +} diff --git a/L1Trigger/DTTriggerPhase2/src/MuonPathFitter.cc b/L1Trigger/DTTriggerPhase2/src/MuonPathFitter.cc new file mode 100644 index 0000000000000..57dc963b79995 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MuonPathFitter.cc @@ -0,0 +1,399 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathFitter.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; + +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MuonPathFitter::MuonPathFitter(const ParameterSet &pset, + edm::ConsumesCollector &iC, + std::shared_ptr &globalcoordsobtainer) + : MuonPathAnalyzer(pset, iC), debug_(pset.getUntrackedParameter("debug")) { + if (debug_) + LogDebug("MuonPathFitter") << "MuonPathAnalyzer: constructor"; + + //shift phi + int rawId; + shift_filename_ = pset.getParameter("shift_filename"); + std::ifstream ifin3(shift_filename_.fullPath()); + double shift; + if (ifin3.fail()) { + throw cms::Exception("Missing Input File") + << "MuonPathFitter::MuonPathFitter() - Cannot find " << shift_filename_.fullPath(); + } + while (ifin3.good()) { + ifin3 >> rawId >> shift; + shiftinfo_[rawId] = shift; + } + + int wh, st, se, maxdrift; + maxdrift_filename_ = pset.getParameter("maxdrift_filename"); + std::ifstream ifind(maxdrift_filename_.fullPath()); + if (ifind.fail()) { + throw cms::Exception("Missing Input File") + << "MPSLFilter::MPSLFilter() - Cannot find " << maxdrift_filename_.fullPath(); + } + while (ifind.good()) { + ifind >> wh >> st >> se >> maxdrift; + maxdriftinfo_[wh][st][se] = maxdrift; + } + + dtGeomH = iC.esConsumes(); + globalcoordsobtainer_ = globalcoordsobtainer; +} + +MuonPathFitter::~MuonPathFitter() { + if (debug_) + LogDebug("MuonPathFitter") << "MuonPathAnalyzer: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +fit_common_out_t MuonPathFitter::fit(fit_common_in_t fit_common_in, + int XI_WIDTH, + int COEFF_WIDTH_T0, + int COEFF_WIDTH_POSITION, + int COEFF_WIDTH_SLOPE, + int PRECISSION_T0, + int PRECISSION_POSITION, + int PRECISSION_SLOPE, + int PROD_RESIZE_T0, + int PROD_RESIZE_POSITION, + int PROD_RESIZE_SLOPE, + int MAX_DRIFT_TDC, + int sl) { + const int PARTIALS_PRECISSION = 4; + const int PARTIALS_SHR_T0 = PRECISSION_T0 - PARTIALS_PRECISSION; + const int PARTIALS_SHR_POSITION = PRECISSION_POSITION - PARTIALS_PRECISSION; + const int PARTIALS_SHR_SLOPE = PRECISSION_SLOPE - PARTIALS_PRECISSION; + const int PARTIALS_WIDTH_T0 = PROD_RESIZE_T0 - PARTIALS_SHR_T0; + const int PARTIALS_WIDTH_POSITION = PROD_RESIZE_POSITION - PARTIALS_SHR_POSITION; + const int PARTIALS_WIDTH_SLOPE = PROD_RESIZE_SLOPE - PARTIALS_SHR_SLOPE; + + const int WIDTH_TO_PREC = 11 + PARTIALS_PRECISSION; + const int WIDTH_SLOPE_PREC = 14 + PARTIALS_PRECISSION; + const int WIDTH_POSITION_PREC = WIDTH_SLOPE_PREC + 1; + + const int SEMICHAMBER_H_PRECISSION = 13 + PARTIALS_PRECISSION; + const float SEMICHAMBER_H_REAL = ((235. / 2.) / (16. * 6.5)) * std::pow(2, SEMICHAMBER_H_PRECISSION); + const int SEMICHAMBER_H = (int)SEMICHAMBER_H_REAL; // signed(SEMICHAMBER_H_WIDTH-1 downto 0) + + const int SEMICHAMBER_RES_SHR = SEMICHAMBER_H_PRECISSION; + + const int LYRANDAHALF_RES_SHR = 4; + + const int CHI2_CALC_RES_BITS = 7; + + /******************************* + clock cycle 1 + *******************************/ + std::vector normalized_times; + std::vector normalized_wirepos; + + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + // normalized times + // this should be resized to an unsigned of 10 bits (max drift time ~508 TDC counts, using 9+1 to include tolerance) + // leaving it as an integer for now + // we are obtaining the difference as the difference in BX + the LS bits from the hit time + + if (fit_common_in.hits_valid[i] == 1) { + int dif_bx = (fit_common_in.hits[i].ti >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME)) - fit_common_in.coarse_bctr; + + int tmp_norm_time = (dif_bx << (WIDTH_FULL_TIME - WIDTH_COARSED_TIME)) + + (fit_common_in.hits[i].ti % (int)std::pow(2, WIDTH_FULL_TIME - WIDTH_COARSED_TIME)); + // resize test + // this has implications in the FW (reducing number of bits). + // we keep here the int as it is, but we do the same check done in the fw + std::vector tmp_dif_bx_vector; + vhdl_int_to_unsigned(dif_bx, tmp_dif_bx_vector); + vhdl_resize_unsigned(tmp_dif_bx_vector, 12); + if (!vhdl_resize_unsigned_ok(tmp_dif_bx_vector, WIDTH_DIFBX)) + return fit_common_out_t(); + + normalized_times.push_back(tmp_norm_time); + int tmp_wirepos = fit_common_in.hits[i].wp - (fit_common_in.coarse_wirepos << WIREPOS_NORM_LSB_IGNORED); + // resize test + std::vector tmp_wirepos_vector; + vhdl_int_to_signed(tmp_wirepos, tmp_wirepos_vector); + vhdl_resize_signed(tmp_wirepos_vector, WIREPOS_WIDTH); + if (!vhdl_resize_signed_ok(tmp_wirepos_vector, XI_WIDTH)) + return fit_common_out_t(); + + normalized_wirepos.push_back(tmp_wirepos); + } else { // dummy hit + normalized_times.push_back(-1); + normalized_wirepos.push_back(-1); + } + } + + /******************************* + clock cycle 2 + *******************************/ + + std::vector xi_arr; + // min and max times are computed throught several clk cycles in the fw, + // here we compute it at once + int min_hit_time = 999999, max_hit_time = 0; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 1) { + // calculate xi array + auto tmp_xi_incr = normalized_wirepos[i]; + tmp_xi_incr += (-1 + 2 * fit_common_in.lateralities[i]) * normalized_times[i]; + + // resize test + std::vector tmp_xi_incr_vector; + vhdl_int_to_signed(tmp_xi_incr, tmp_xi_incr_vector); + vhdl_resize_signed(tmp_xi_incr_vector, XI_WIDTH + 1); + if (!vhdl_resize_signed_ok(tmp_xi_incr_vector, XI_WIDTH)) + return fit_common_out_t(); + xi_arr.push_back(tmp_xi_incr); + + // calculate min and max times + if (normalized_times[i] < min_hit_time) { + min_hit_time = normalized_times[i]; + } + if (normalized_times[i] > max_hit_time) { + max_hit_time = normalized_times[i]; + } + } else { + xi_arr.push_back(-1); + } + } + + /******************************* + clock cycle 3 + *******************************/ + + std::vector products_t0; + std::vector products_position; + std::vector products_slope; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 0) { + products_t0.push_back(-1); + products_position.push_back(-1); + products_slope.push_back(-1); + } else { + products_t0.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.t0[i])); + products_position.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.position[i])); + products_slope.push_back(xi_arr[i] * vhdl_signed_to_int(fit_common_in.coeffs.slope[i])); + } + } + + /******************************* + clock cycle 4 + *******************************/ + // Do the 8 element sums + int t0_prec = 0, position_prec = 0, slope_prec = 0; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 0) { + continue; + } else { + t0_prec += products_t0[i] >> PARTIALS_SHR_T0; + position_prec += products_position[i] >> PARTIALS_SHR_POSITION; + slope_prec += products_slope[i] >> PARTIALS_SHR_SLOPE; + } + } + + /******************************* + clock cycle 5 + *******************************/ + // Do resize tests for the computed sums with full precision + std::vector t0_prec_vector, position_prec_vector, slope_prec_vector; + vhdl_int_to_signed(t0_prec, t0_prec_vector); + + vhdl_resize_signed(t0_prec_vector, PARTIALS_WIDTH_T0); + if (!vhdl_resize_signed_ok(t0_prec_vector, WIDTH_TO_PREC)) + return fit_common_out_t(); + + vhdl_int_to_signed(position_prec, position_prec_vector); + vhdl_resize_signed(position_prec_vector, PARTIALS_WIDTH_POSITION); + if (!vhdl_resize_signed_ok(position_prec_vector, WIDTH_POSITION_PREC)) + return fit_common_out_t(); + + vhdl_int_to_signed(slope_prec, slope_prec_vector); + vhdl_resize_signed(slope_prec_vector, PARTIALS_WIDTH_SLOPE); + if (!vhdl_resize_signed_ok(slope_prec_vector, WIDTH_SLOPE_PREC)) + return fit_common_out_t(); + + /******************************* + clock cycle 6 + *******************************/ + // Round the fitting parameters to the final resolution; + // in vhdl something more sofisticated is done, here we do a float division, round + // and cast again to integer + + int norm_t0 = ((t0_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1; + int norm_position = ((position_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1; + int norm_slope = ((slope_prec >> (PARTIALS_PRECISSION - 1)) + 1) >> 1; + + // Calculate the (-xi) + pos (+/-) t0, which only is lacking the slope term to become the residuals + std::vector res_partials_arr; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 0) { + res_partials_arr.push_back(-1); + } else { + int tmp_position_prec = position_prec - (xi_arr[i] << PARTIALS_PRECISSION); + // rounding + tmp_position_prec += std::pow(2, PARTIALS_PRECISSION - 1); + + tmp_position_prec += (-1 + 2 * fit_common_in.lateralities[i]) * t0_prec; + res_partials_arr.push_back(tmp_position_prec); + } + } + + // calculate the { slope x semichamber, slope x 1.5 layers, slope x 0.5 layers } + // these 3 values are later combined with different signs to get the slope part + // of the residual for each of the layers. + int slope_x_halfchamb = (((long int)slope_prec * (long int)SEMICHAMBER_H)) >> SEMICHAMBER_RES_SHR; + if (sl == 2) + slope_x_halfchamb = 0; + int slope_x_3semicells = (slope_prec * 3) >> LYRANDAHALF_RES_SHR; + int slope_x_1semicell = (slope_prec * 1) >> LYRANDAHALF_RES_SHR; + + /******************************* + clock cycle 7 + *******************************/ + // Complete the residuals calculation by constructing the slope term (1/2) + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 1) { + if (i % 4 == 0) + res_partials_arr[i] -= slope_x_3semicells; + else if (i % 4 == 1) + res_partials_arr[i] -= slope_x_1semicell; + else if (i % 4 == 2) + res_partials_arr[i] += slope_x_1semicell; + else + res_partials_arr[i] += slope_x_3semicells; + } + } + + /******************************* + clock cycle 8 + *******************************/ + // Complete the residuals calculation by constructing the slope term (2/2) + std::vector residuals, position_prec_arr; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 0) { + residuals.push_back(-1); + position_prec_arr.push_back(-1); + } else { + int tmp_position_prec = res_partials_arr[i]; + tmp_position_prec += (-1 + 2 * (int)(i >= NUM_LAYERS)) * slope_x_halfchamb; + position_prec_arr.push_back(tmp_position_prec); + residuals.push_back(abs(tmp_position_prec >> PARTIALS_PRECISSION)); + } + } + + // minimum and maximum fit t0 + int min_t0 = max_hit_time - MAX_DRIFT_TDC - T0_CUT_TOLERANCE; + int max_t0 = min_hit_time + T0_CUT_TOLERANCE; + + /******************************* + clock cycle 9 + *******************************/ + // Prepare addition of coarse_offset to T0 (T0 de-normalization) + int t0_fine = norm_t0 & (int)(std::pow(2, 5) - 1); + int t0_bx_sign = ((int)(norm_t0 < 0)) * 1; + int t0_bx_abs = abs(norm_t0 >> 5); + + // De-normalize Position and slope + int position = (fit_common_in.coarse_wirepos << WIREPOS_NORM_LSB_IGNORED) + norm_position; + int slope = norm_slope; + + // Apply T0 cuts + if (norm_t0 < min_t0) + return fit_common_out_t(); + if (norm_t0 > max_t0) + return fit_common_out_t(); + + // square the residuals + std::vector squared_residuals; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 0) { + squared_residuals.push_back(-1); + } else { + squared_residuals.push_back(residuals[i] * residuals[i]); + } + } + + // check for residuals overflow + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 1) { + std::vector tmp_vector; + int tmp_position_prec = (position_prec_arr[i] >> PARTIALS_PRECISSION); + vhdl_int_to_signed(tmp_position_prec, tmp_vector); + vhdl_resize_signed(tmp_vector, WIDTH_POSITION_PREC); + if (!vhdl_resize_signed_ok(tmp_vector, CHI2_CALC_RES_BITS + 1)) + return fit_common_out_t(); + // Commented for now, maybe later we need to do something here + // if ((tmp_position_prec / (int) std::pow(2, CHI2_CALC_RES_BITS)) > 0) + // return fit_common_out_t(); + } + } + + /******************************* + clock cycle 10, 11, 12 + *******************************/ + int t0 = t0_fine; + t0 += (fit_common_in.coarse_bctr - (-1 + 2 * t0_bx_sign) * t0_bx_abs) * (int)std::pow(2, 5); + + int chi2 = 0; + for (int i = 0; i < 2 * NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[i] == 1) { + chi2 += squared_residuals[i]; + } + } + + // Impose the thresholds + + if (chi2 / 16 >= (int)round(chi2Th_ * (std::pow((float)MAX_DRIFT_TDC / ((float)CELL_SEMILENGTH / 10.), 2)) / 16.)) + return fit_common_out_t(); + + fit_common_out_t fit_common_out; + fit_common_out.position = position; + fit_common_out.slope = slope; + fit_common_out.t0 = t0; + fit_common_out.chi2 = chi2; + fit_common_out.valid_fit = 1; + + return fit_common_out; +} + +coeffs_t MuonPathFitter::RomDataConvert(std::vector slv, + short COEFF_WIDTH_T0, + short COEFF_WIDTH_POSITION, + short COEFF_WIDTH_SLOPE, + short LOLY, + short HILY) { + coeffs_t res; + int ctr = 0; + for (int i = LOLY; i <= HILY; i++) { + res.t0[i] = vhdl_slice(slv, COEFF_WIDTH_T0 + ctr - 1, ctr); + vhdl_resize_unsigned(res.t0[i], GENERIC_COEFF_WIDTH); + res.t0[i] = vhdl_slice(res.t0[i], COEFF_WIDTH_T0 - 1, 0); + ctr += COEFF_WIDTH_T0; + } + for (int i = LOLY; i <= HILY; i++) { + res.position[i] = vhdl_slice(slv, COEFF_WIDTH_POSITION + ctr - 1, ctr); + vhdl_resize_unsigned(res.position[i], GENERIC_COEFF_WIDTH); + res.position[i] = vhdl_slice(res.position[i], COEFF_WIDTH_POSITION - 1, 0); + ctr += COEFF_WIDTH_POSITION; + } + for (int i = LOLY; i <= HILY; i++) { + res.slope[i] = vhdl_slice(slv, COEFF_WIDTH_SLOPE + ctr - 1, ctr); + vhdl_resize_unsigned(res.slope[i], GENERIC_COEFF_WIDTH); + res.slope[i] = vhdl_slice(res.slope[i], COEFF_WIDTH_SLOPE - 1, 0); + ctr += COEFF_WIDTH_SLOPE; + } + return res; +} diff --git a/L1Trigger/DTTriggerPhase2/src/MuonPathSLFitter.cc b/L1Trigger/DTTriggerPhase2/src/MuonPathSLFitter.cc new file mode 100644 index 0000000000000..fda1c83d9cdaa --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/MuonPathSLFitter.cc @@ -0,0 +1,509 @@ +#include "L1Trigger/DTTriggerPhase2/interface/MuonPathSLFitter.h" +#include +#include + +using namespace edm; +using namespace std; +using namespace cmsdt; +// ============================================================================ +// Constructors and destructor +// ============================================================================ +MuonPathSLFitter::MuonPathSLFitter(const ParameterSet &pset, + edm::ConsumesCollector &iC, + std::shared_ptr &globalcoordsobtainer) + : MuonPathFitter(pset, iC, globalcoordsobtainer) { + if (debug_) + LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: constructor"; + + //shift theta + int rawId; + double shift; + shift_theta_filename_ = pset.getParameter("shift_theta_filename"); + std::ifstream ifin4(shift_theta_filename_.fullPath()); + if (ifin4.fail()) { + throw cms::Exception("Missing Input File") + << "MuonPathSLFitter::MuonPathSLFitter() - Cannot find " << shift_theta_filename_.fullPath(); + } + + while (ifin4.good()) { + ifin4 >> rawId >> shift; + shiftthetainfo_[rawId] = shift; + } + + // LUTs + sl1_filename_ = pset.getParameter("lut_sl1"); + sl2_filename_ = pset.getParameter("lut_sl2"); + sl3_filename_ = pset.getParameter("lut_sl3"); + + fillLuts(); + + setChi2Th(pset.getParameter("chi2Th")); + setTanPhiTh(pset.getParameter("tanPhiTh")); +} + +MuonPathSLFitter::~MuonPathSLFitter() { + if (debug_) + LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void MuonPathSLFitter::initialise(const edm::EventSetup &iEventSetup) { + if (debug_) + LogDebug("MuonPathSLFitter") << "MuonPathSLFitter::initialiase"; + + auto geom = iEventSetup.getHandle(dtGeomH); + dtGeo_ = &(*geom); +} + +void MuonPathSLFitter::run(edm::Event &iEvent, + const edm::EventSetup &iEventSetup, + MuonPathPtrs &muonpaths, + std::vector &lateralities, + std::vector &metaPrimitives) { + if (debug_) + LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: run"; + + // fit per SL (need to allow for multiple outputs for a single mpath) + // for (auto &muonpath : muonpaths) { + if (!muonpaths.empty()) { + auto muonpath = muonpaths[0]; + int rawId = muonpath->primitive(0)->cameraId(); + if (muonpath->primitive(0)->cameraId() == -1) { + rawId = muonpath->primitive(1)->cameraId(); + } + const DTLayerId dtLId(rawId); + max_drift_tdc = maxdriftinfo_[dtLId.wheel() + 2][dtLId.station() - 1][dtLId.sector() - 1]; + } + + for (size_t i = 0; i < muonpaths.size(); i++) { + auto muonpath = muonpaths[i]; + auto lats = lateralities[i]; + analyze(muonpath, lats, metaPrimitives); + } +} + +void MuonPathSLFitter::finish() { + if (debug_) + LogDebug("MuonPathSLFitter") << "MuonPathSLFitter: finish"; +}; + +//------------------------------------------------------------------ +//--- Metodos privados +//------------------------------------------------------------------ + +void MuonPathSLFitter::analyze(MuonPathPtr &inMPath, + lat_vector lat_combs, + std::vector &metaPrimitives) { + auto sl = inMPath->primitive(0)->superLayerId(); // 0, 1, 2 + + int selected_lay = 1; + if (inMPath->primitive(0)->tdcTimeStamp() != -1) + selected_lay = 0; + + int dumLayId = inMPath->primitive(selected_lay)->cameraId(); + auto dtDumlayerId = DTLayerId(dumLayId); + DTSuperLayerId MuonPathSLId(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), sl + 1); + + DTChamberId ChId(MuonPathSLId.wheel(), MuonPathSLId.station(), MuonPathSLId.sector()); + + DTSuperLayerId MuonPathSL1Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 1); + DTSuperLayerId MuonPathSL2Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 2); + DTSuperLayerId MuonPathSL3Id(dtDumlayerId.wheel(), dtDumlayerId.station(), dtDumlayerId.sector(), 3); + DTWireId wireIdSL1(MuonPathSL1Id, 2, 1); + DTWireId wireIdSL2(MuonPathSL2Id, 2, 1); + DTWireId wireIdSL3(MuonPathSL3Id, 2, 1); + auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()]; + + fit_common_in_t fit_common_in; + + // 8-element vectors, for the 8 layers. As here we are fitting one SL only, we leave the other SL values as dummy ones + fit_common_in.hits = {}; + fit_common_in.hits_valid = {}; + + int quality = 3; + if (inMPath->missingLayer() != -1) + quality = 1; + + int minISL = 1; + int maxISL = 3; + if (sl < 1) { + minISL = 0; + maxISL = 2; + } + + for (int isl = minISL; isl < maxISL; isl++) { + for (int i = 0; i < NUM_LAYERS; i++) { + if (isl == sl && inMPath->missingLayer() != i) { + // Include both valid and non-valid hits. Non-valid values can be whatever, leaving all as -1 to make debugging easier. + auto ti = inMPath->primitive(i)->tdcTimeStamp(); + if (ti != -1) + ti = (int)round(((float)TIME_TO_TDC_COUNTS / (float)LHC_CLK_FREQ) * ti); + auto wi = inMPath->primitive(i)->channelId(); + auto ly = inMPath->primitive(i)->layerId(); + // int layId = inMPath->primitive(i)->cameraId(); + // auto dtlayerId = DTLayerId(layId); + // auto wireId = DTWireId(dtlayerId, wi + 1); // wire start from 1, mixer groups them starting from 0 + // int rawId = wireId.rawId(); + // wp in tdc counts (still in floating point) + int wp_semicells = (wi - SL1_CELLS_OFFSET) * 2 + 1; + if (ly % 2 == 1) + wp_semicells -= 1; + if (isl == 2) + wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH); + + float wp_tdc = wp_semicells * max_drift_tdc; + int wp = (int)((long int)(round(wp_tdc * std::pow(2, WIREPOS_WIDTH))) / (int)std::pow(2, WIREPOS_WIDTH)); + fit_common_in.hits.push_back({ti, wi, ly, wp}); + // fill valids as well + if (inMPath->missingLayer() == i) + fit_common_in.hits_valid.push_back(0); + else + fit_common_in.hits_valid.push_back(1); + } else { + fit_common_in.hits.push_back({-1, -1, -1, -1}); + fit_common_in.hits_valid.push_back(0); + } + } + } + + int smallest_time = 999999, tmp_coarse_wirepos_1 = -1, tmp_coarse_wirepos_3 = -1; + // coarse_bctr is the 12 MSB of the smallest tdc + for (int isl = 0; isl < 3; isl++) { + if (isl != sl) + continue; + int myisl = isl < 2 ? 0 : 1; + for (size_t i = 0; i < NUM_LAYERS; i++) { + if (fit_common_in.hits_valid[NUM_LAYERS * myisl + i] == 0) + continue; + else if (fit_common_in.hits[NUM_LAYERS * myisl + i].ti < smallest_time) + smallest_time = fit_common_in.hits[NUM_LAYERS * myisl + i].ti; + } + if (fit_common_in.hits_valid[NUM_LAYERS * myisl + 0] == 1) + tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * myisl + 0].wp; + else + tmp_coarse_wirepos_1 = fit_common_in.hits[NUM_LAYERS * myisl + 1].wp; + if (fit_common_in.hits_valid[NUM_LAYERS * myisl + 3] == 1) + tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * myisl + 3].wp; + else + tmp_coarse_wirepos_3 = fit_common_in.hits[NUM_LAYERS * myisl + 2].wp; + + tmp_coarse_wirepos_1 = tmp_coarse_wirepos_1 >> WIREPOS_NORM_LSB_IGNORED; + tmp_coarse_wirepos_3 = tmp_coarse_wirepos_3 >> WIREPOS_NORM_LSB_IGNORED; + } + fit_common_in.coarse_bctr = smallest_time >> (WIDTH_FULL_TIME - WIDTH_COARSED_TIME); + fit_common_in.coarse_wirepos = (tmp_coarse_wirepos_1 + tmp_coarse_wirepos_3) >> 1; + + for (auto &lat_comb : lat_combs) { + if (lat_comb[0] == 0 && lat_comb[1] == 0 && lat_comb[2] == 0 && lat_comb[3] == 0) + continue; + fit_common_in.lateralities.clear(); + + auto rom_addr = get_rom_addr(inMPath, lat_comb); + coeffs_t coeffs; + if (sl == 0) { + coeffs = + RomDataConvert(lut_sl1[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL_POSITION, COEFF_WIDTH_SL_SLOPE, 0, 3); + } else if (sl == 1) { + coeffs = + RomDataConvert(lut_sl2[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL2_POSITION, COEFF_WIDTH_SL_SLOPE, 0, 3); + } else { + coeffs = + RomDataConvert(lut_sl3[rom_addr], COEFF_WIDTH_SL_T0, COEFF_WIDTH_SL_POSITION, COEFF_WIDTH_SL_SLOPE, 4, 7); + } + // Filling lateralities + int minISL = 1; + int maxISL = 3; + if (sl < 1) { + minISL = 0; + maxISL = 2; + } + + for (int isl = minISL; isl < maxISL; isl++) { + for (size_t i = 0; i < NUM_LAYERS; i++) { + if (isl == sl) { + fit_common_in.lateralities.push_back(lat_comb[i]); + } else + fit_common_in.lateralities.push_back(-1); + } + } + fit_common_in.coeffs = coeffs; + + auto fit_common_out = fit(fit_common_in, + XI_SL_WIDTH, + COEFF_WIDTH_SL_T0, + sl == 1 ? COEFF_WIDTH_SL2_POSITION : COEFF_WIDTH_SL_POSITION, + COEFF_WIDTH_SL_SLOPE, + PRECISSION_SL_T0, + PRECISSION_SL_POSITION, + PRECISSION_SL_SLOPE, + PROD_RESIZE_SL_T0, + sl == 1 ? PROD_RESIZE_SL2_POSITION : PROD_RESIZE_SL_POSITION, + PROD_RESIZE_SL_SLOPE, + max_drift_tdc, + sl + 1); + + if (fit_common_out.valid_fit == 1) { + float t0_f = ((float)fit_common_out.t0) * (float)LHC_CLK_FREQ / (float)TIME_TO_TDC_COUNTS; + + float slope_f = -fit_common_out.slope * ((float)CELL_SEMILENGTH / max_drift_tdc) * (1) / (CELL_SEMIHEIGHT * 16.); + if (sl != 1 && std::abs(slope_f) > tanPhiTh_) + continue; + + DTWireId wireId(MuonPathSLId, 2, 1); + float pos_ch_f = (float)(fit_common_out.position) * ((float)CELL_SEMILENGTH / (float)max_drift_tdc) / 10; + pos_ch_f += (SL1_CELLS_OFFSET * CELL_LENGTH) / 10.; + if (sl != 1) + pos_ch_f += shiftinfo_[wireIdSL1.rawId()]; + else if (sl == 1) + pos_ch_f += shiftthetainfo_[wireIdSL2.rawId()]; + + float pos_sl_f = pos_ch_f - (sl - 1) * slope_f * VERT_PHI1_PHI3 / 2; + float chi2_f = fit_common_out.chi2 * std::pow(((float)CELL_SEMILENGTH / (float)max_drift_tdc), 2) / 100; + + // obtention of global coordinates using luts + int pos = (int)(10 * (pos_sl_f - shiftinfo_[wireId.rawId()]) * INCREASED_RES_POS_POW); + int slope = (int)(-slope_f * INCREASED_RES_SLOPE_POW); + + auto global_coords = globalcoordsobtainer_->get_global_coordinates(ChId.rawId(), sl + 1, pos, slope); + float phi = global_coords[0]; + float phiB = global_coords[1]; + + // obtention of global coordinates using cmssw geometry + double z = 0; + if (ChId.station() == 3 or ChId.station() == 4) { + z = Z_SHIFT_MB4; + } + GlobalPoint jm_x_cmssw_global = dtGeo_->chamber(ChId)->toGlobal(LocalPoint(pos_sl_f, 0., z)); + int thisec = ChId.sector(); + if (thisec == 13) + thisec = 4; + if (thisec == 14) + thisec = 10; + float phi_cmssw = jm_x_cmssw_global.phi() - PHI_CONV * (thisec - 1); + float psi = atan(slope_f); + float phiB_cmssw = hasPosRF(ChId.wheel(), ChId.sector()) ? psi - phi_cmssw : -psi - phi_cmssw; + if (sl == 0) + metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(), + t0_f, + (double)(fit_common_out.position), + (double)fit_common_out.slope, + phi, + phiB, + phi_cmssw, + phiB_cmssw, + chi2_f, + quality, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + lat_comb[0], + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + lat_comb[1], + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + lat_comb[2], + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + lat_comb[3], + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1})); + else if (sl == 2) + metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(), + t0_f, + (double)(fit_common_out.position), + (double)fit_common_out.slope, + phi, + phiB, + phi_cmssw, + phiB_cmssw, + chi2_f, + quality, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + lat_comb[0], + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + lat_comb[1], + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + lat_comb[2], + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + lat_comb[3], + -1})); + else if (sl == 1) { + // fw-like calculation + DTLayerId SL2_layer2Id(MuonPathSLId, 2); + double z_shift = shiftthetainfo_[SL2_layer2Id.rawId()]; + double pos_cm = + pos / 10 / INCREASED_RES_POS_POW; // fixed to have z_shift and the position in the same units (MC) + double jm_y = hasPosRF(MuonPathSLId.wheel(), MuonPathSLId.sector()) ? z_shift - pos_cm : z_shift + pos_cm; + + phi = jm_y; + + // Fixed sign of k (MC) + double k_fromfw = hasPosRF(MuonPathSLId.wheel(), MuonPathSLId.sector()) ? slope_f : -slope_f; + phiB = k_fromfw; + + // cmssw-like calculation + LocalPoint wire1_in_layer(dtGeo_->layer(SL2_layer2Id)->specificTopology().wirePosition(1), 0, -0.65); + GlobalPoint wire1_in_global = dtGeo_->layer(SL2_layer2Id)->toGlobal(wire1_in_layer); + LocalPoint wire1_in_sl = dtGeo_->superLayer(MuonPathSLId)->toLocal(wire1_in_global); + double x_shift = wire1_in_sl.x(); + jm_y = (dtGeo_->superLayer(MuonPathSLId) + ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0))) + .z(); + phi_cmssw = jm_y; + + double x_cmssw = (dtGeo_->superLayer(MuonPathSLId) + ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0))) + .x(); + double y_cmssw = (dtGeo_->superLayer(MuonPathSLId) + ->toGlobal(LocalPoint(double(pos) / (10 * pow(2, INCREASED_RES_POS)) + x_shift, 0., 0))) + .y(); + double r_cmssw = sqrt(x_cmssw * x_cmssw + y_cmssw * y_cmssw); + double k_cmssw = jm_y / r_cmssw; + phiB_cmssw = k_cmssw; + + metaPrimitives.emplace_back(metaPrimitive({MuonPathSLId.rawId(), + t0_f, + (double)(fit_common_out.position), + (double)fit_common_out.slope, + phi, + phiB, + phi_cmssw, + phiB_cmssw, + chi2_f, + quality, + inMPath->primitive(0)->channelId(), + inMPath->primitive(0)->tdcTimeStamp(), + lat_comb[0], + inMPath->primitive(1)->channelId(), + inMPath->primitive(1)->tdcTimeStamp(), + lat_comb[1], + inMPath->primitive(2)->channelId(), + inMPath->primitive(2)->tdcTimeStamp(), + lat_comb[2], + inMPath->primitive(3)->channelId(), + inMPath->primitive(3)->tdcTimeStamp(), + lat_comb[3], + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1, + -1})); + } + } // (fit_common_out.valid_fit == 1) + } // loop in lat_combs + return; +} + +void MuonPathSLFitter::fillLuts() { + std::ifstream ifinsl1(sl1_filename_.fullPath()); + std::string line; + while (ifinsl1.good()) { + ifinsl1 >> line; + + std::vector myNumbers; + for (size_t i = 0; i < line.size(); i++) { + // This converts the char into an int and pushes it into vec + myNumbers.push_back(line[i] - '0'); // The digits will be in the same order as before + } + std::reverse(myNumbers.begin(), myNumbers.end()); + lut_sl1.push_back(myNumbers); + } + + std::ifstream ifinsl2(sl2_filename_.fullPath()); + while (ifinsl2.good()) { + ifinsl2 >> line; + + std::vector myNumbers; + for (size_t i = 0; i < line.size(); i++) { + // This converts the char into an int and pushes it into vec + myNumbers.push_back(line[i] - '0'); // The digits will be in the same order as before + } + std::reverse(myNumbers.begin(), myNumbers.end()); + lut_sl2.push_back(myNumbers); + } + + std::ifstream ifinsl3(sl3_filename_.fullPath()); + while (ifinsl3.good()) { + ifinsl3 >> line; + + std::vector myNumbers; + for (size_t i = 0; i < line.size(); i++) { + // This converts the char into an int and pushes it into vec + myNumbers.push_back(line[i] - '0'); // The digits will be in the same order as before + } + std::reverse(myNumbers.begin(), myNumbers.end()); + lut_sl3.push_back(myNumbers); + } + + return; +} + +int MuonPathSLFitter::get_rom_addr(MuonPathPtr &inMPath, latcomb lats) { + std::vector rom_addr; + auto missing_layer = inMPath->missingLayer(); + if (missing_layer == -1) { + rom_addr.push_back(1); + rom_addr.push_back(0); + } else { + if (missing_layer == 0) { + rom_addr.push_back(0); + rom_addr.push_back(0); + } else if (missing_layer == 1) { + rom_addr.push_back(0); + rom_addr.push_back(1); + } else if (missing_layer == 2) { + rom_addr.push_back(1); + rom_addr.push_back(0); + } else { // missing_layer == 3 + rom_addr.push_back(1); + rom_addr.push_back(1); + } + } + for (size_t ilat = 0; ilat < lats.size(); ilat++) { + if ((int)ilat == missing_layer) // only applies to 3-hit, as in 4-hit missL=-1 + continue; + auto lat = lats[ilat]; + if (lat == -1) + lat = 0; + rom_addr.push_back(lat); + } + std::reverse(rom_addr.begin(), rom_addr.end()); + return vhdl_unsigned_to_int(rom_addr); +} diff --git a/L1Trigger/DTTriggerPhase2/src/TrapezoidalGrouping.cc b/L1Trigger/DTTriggerPhase2/src/TrapezoidalGrouping.cc new file mode 100644 index 0000000000000..fee10c32095d0 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/TrapezoidalGrouping.cc @@ -0,0 +1,229 @@ +#include + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h" + +using namespace edm; +using namespace std; +using namespace cmsdt; +using namespace dtamgrouping; +// ============================================================================ +// Constructors and destructor +// ============================================================================ +TrapezoidalGrouping::TrapezoidalGrouping(const ParameterSet &pset, edm::ConsumesCollector &iC) + : MotherGrouping(pset, iC), debug_(pset.getUntrackedParameter("debug")), currentBaseChannel_(-1) { + // Obtention of parameters + if (debug_) + LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping: constructor"; + + // Initialisation of channelIn array + chInDummy_.push_back(DTPrimitive()); + for (int lay = 0; lay < NUM_LAYERS; lay++) { + for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) { + channelIn_[lay][ch] = {chInDummy_}; + channelIn_[lay][ch].clear(); + } + } +} + +TrapezoidalGrouping::~TrapezoidalGrouping() { + if (debug_) + LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping: destructor"; +} + +// ============================================================================ +// Main methods (initialise, run, finish) +// ============================================================================ +void TrapezoidalGrouping::initialise(const edm::EventSetup &iEventSetup) { + if (debug_) + LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping::initialiase"; +} + +void TrapezoidalGrouping::run(Event &iEvent, + const EventSetup &iEventSetup, + const DTDigiCollection &digis, + MuonPathPtrs &mpaths) { + // This function returns the analyzable mpath collection back to the the main function + // so it can be fitted. This is in fact doing the so-called grouping. + for (int supLayer = 0; supLayer < NUM_SUPERLAYERS; supLayer++) { // for each SL: + if (debug_) + LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping::run Reading SL" << supLayer; + setInChannels(&digis, supLayer); + + for (auto &hit : all_hits) { + int layer_to_pivot = hit.layerId(); + int channel_to_pivot = hit.channelId(); + DTPrimitives hits_in_trapezoid; + std::vector hit_mpaths; + std::vector hit_tasks; + for (size_t itask = 0; itask < task_list.size(); itask++) { + // when pivoting over an internal layer, there are two cases + // where the second layer is duplicated + // 12 (0, 5) <-> 14 (0, 7) + // 15 (1, 6) <-> 17 (1, 8) + // we leave it hard-coded here, could be moved somewhere else + if (layer_to_pivot == 1 || layer_to_pivot == 2) { + if (itask == 14 || itask == 17) + continue; + } + + auto task = task_list[itask]; + + std::vector task_mpaths; + std::stack> mpath_cells_per_task; + mpath_cells_per_task.push(std::make_pair(DTPrimitives({hit}), 0)); + + while (!mpath_cells_per_task.empty()) { + auto mpath_cells = std::move(mpath_cells_per_task.top()); + std::vector tmp_mpaths = {mpath_cells.first}; + auto task_index = mpath_cells.second; + auto cell = task[task_index]; + auto vertical_shift = trapezoid_vertical_mapping[layer_to_pivot][cell]; + auto horizontal_shift = trapezoid_horizontal_mapping[layer_to_pivot][cell]; + if (channel_to_pivot + horizontal_shift >= 0 && channel_to_pivot + horizontal_shift < NUM_CH_PER_LAYER) { + tmp_mpaths = group_hits(hit, + tmp_mpaths, + channelIn_[layer_to_pivot + vertical_shift][channel_to_pivot + horizontal_shift], + hits_in_trapezoid); + } + mpath_cells_per_task.pop(); + for (const auto &tmp_mpath : tmp_mpaths) { + mpath_cells_per_task.push(std::make_pair(tmp_mpath, task_index + 1)); + } + while (!mpath_cells_per_task.empty()) { + if (mpath_cells_per_task.top().second == (int)task.size()) { + task_mpaths.push_back(mpath_cells_per_task.top().first); + mpath_cells_per_task.pop(); + } else + break; + } + } + for (auto &task_mpath : task_mpaths) { + hit_mpaths.push_back(task_mpath); + hit_tasks.push_back(itask); + } + } + if (hits_in_trapezoid.size() <= PATHFINDER_INPUT_HITS_LIMIT) { + for (size_t ipath = 0; ipath < hit_mpaths.size(); ipath++) { + auto ptrPrimitive = hit_mpaths[ipath]; + auto itask = hit_tasks[ipath]; + + // In any case, if we have less than 3 hits, we don't output the mpath + if (ptrPrimitive.size() < 3) + continue; + + // check if the task has a missing layer associated + // if it does, we add a dummy hit in the missing layer + // if it does not, we check that we actually have 4 present hits; + // if not, we skip the mpath. + if (MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask] != -1) { + auto dtpAux = DTPrimitive(); + dtpAux.setTDCTimeStamp(-1); + dtpAux.setChannelId(-1); + dtpAux.setLayerId(MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask]); // L=0,1,2,3 + dtpAux.setSuperLayerId(hit.superLayerId()); + dtpAux.setCameraId(-1); + ptrPrimitive.push_back(dtpAux); + } else { // we have no missing hits, it must be a 4-hit TP. + if (ptrPrimitive.size() < 4) + continue; + } + + // sort the hits by layer, so they are included ordered in the MuonPath object + std::stable_sort(ptrPrimitive.begin(), ptrPrimitive.end(), hitLayerSort); + + auto ptrMuonPath = std::make_shared(ptrPrimitive); + ptrMuonPath->setCellHorizontalLayout(CELL_HORIZONTAL_LAYOUTS_PER_TASK[layer_to_pivot][itask]); + ptrMuonPath->setMissingLayer(MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask]); + mpaths.push_back(std::move(ptrMuonPath)); + } + } + } + } + if (debug_) + LogDebug("TrapezoidalGrouping") << "[TrapezoidalGrouping::run] end"; +} + +void TrapezoidalGrouping::finish() { return; }; + +// ============================================================================ +// Other methods +// ============================================================================ + +void TrapezoidalGrouping::setInChannels(const DTDigiCollection *digis, int sl) { + // before setting channels we need to clear + for (int lay = 0; lay < NUM_LAYERS; lay++) { + for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) { + channelIn_[lay][ch].clear(); + } + } + all_hits.clear(); + + // now fill with those primitives that make sense: + for (const auto &dtLayerId_It : *digis) { + const DTLayerId dtLId = dtLayerId_It.first; + + if (dtLId.superlayer() != sl + 1) + continue; //skip digis not in SL... + + for (DTDigiCollection::const_iterator digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; + ++digiIt) { + int layer = dtLId.layer() - 1; + int wire = (*digiIt).wire() - 1; + int digiTIME = (*digiIt).time(); + int digiTIMEPhase2 = digiTIME; + + if (debug_) + LogDebug("TrapezoidalGrouping") << "[TrapezoidalGrouping::setInChannels] SL" << sl << " L" << layer << " : " + << wire << " " << digiTIMEPhase2; + auto dtpAux = DTPrimitive(); + dtpAux.setTDCTimeStamp(digiTIMEPhase2); + dtpAux.setChannelId(wire); + dtpAux.setLayerId(layer); // L=0,1,2,3 + dtpAux.setSuperLayerId(sl); // SL=0,1,2 + dtpAux.setCameraId(dtLId.rawId()); + channelIn_[layer][wire].push_back(dtpAux); + all_hits.push_back(dtpAux); + } + } + + // sort everything by the time of the hits, so it has the same behaviour as the fw + for (int lay = 0; lay < NUM_LAYERS; lay++) { + for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) { + std::stable_sort(channelIn_[lay][ch].begin(), channelIn_[lay][ch].end(), hitTimeSort); + } + } + std::stable_sort(all_hits.begin(), all_hits.end(), hitTimeSort); +} + +std::vector TrapezoidalGrouping::group_hits(DTPrimitive pivot_hit, + std::vector input_paths, + DTPrimitives hits_per_cell, + DTPrimitives &hits_in_trapezoid) { + std::vector output_paths; + for (auto &hit : hits_per_cell) { + int hit_bx = hit.tdcTimeStamp() / LHC_CLK_FREQ; + int pivot_hit_bx = pivot_hit.tdcTimeStamp() / LHC_CLK_FREQ; + if (hitTimeSort(pivot_hit, hit) || (pivot_hit_bx / BX_PER_FRAME) - (hit_bx / BX_PER_FRAME) > MAX_FRAME_DIF) + continue; + // limit the number of hits in the trapezoid to PATHFINDER_INPUT_HITS_LIMIT + if (std::find(hits_in_trapezoid.begin(), hits_in_trapezoid.end(), hit) == hits_in_trapezoid.end()) + hits_in_trapezoid.push_back(hit); + + if (hits_in_trapezoid.size() > PATHFINDER_INPUT_HITS_LIMIT) { + std::vector empty_paths; + return empty_paths; + } + + for (auto &input_path : input_paths) { + auto tmp_path = input_path; + tmp_path.push_back(hit); + output_paths.push_back(tmp_path); + } + } + if (output_paths.empty()) + return input_paths; + else + return output_paths; +} diff --git a/L1Trigger/DTTriggerPhase2/src/lat_code_maker.py b/L1Trigger/DTTriggerPhase2/src/lat_code_maker.py new file mode 100644 index 0000000000000..23a825e65786c --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/lat_code_maker.py @@ -0,0 +1,9 @@ +filename = "basic_lats.txt" + +with open(filename) as f: + lines = f.readlines() + +for iline, line in enumerate(lines): + line = line.strip().split("/") + text = "lat_combinations.push_back({ %s, %s, %s });" % (line[0], line[1], line[2]) + print(text) diff --git a/L1Trigger/DTTriggerPhase2/src/vhdl_functions.cc b/L1Trigger/DTTriggerPhase2/src/vhdl_functions.cc new file mode 100644 index 0000000000000..1343ed5786fb1 --- /dev/null +++ b/L1Trigger/DTTriggerPhase2/src/vhdl_functions.cc @@ -0,0 +1,93 @@ +#include "L1Trigger/DTTriggerPhase2/interface/vhdl_functions.h" + +// "a la vhdl" functions +std::vector vhdl_slice(std::vector v, int upper, int lower) { + int final_value = lower; + if (final_value < 0) + final_value = 0; + + std::vector v1; + for (int i = final_value; i <= upper; i++) { + v1.push_back(v[i]); + } + return v1; +} + +int vhdl_unsigned_to_int(std::vector v) { + int res = 0; + + for (size_t i = 0; i < v.size(); i++) { + res = res + v[i] * std::pow(2, i); + } + return res; +} + +int vhdl_signed_to_int(std::vector v) { + if (v[v.size() - 1] == 0) + return vhdl_unsigned_to_int(v); + else + return -(std::pow(2, v.size()) - vhdl_unsigned_to_int(v)); +} + +void vhdl_int_to_unsigned(int value, std::vector &v) { + if (value == 0) { + v.push_back(0); + } else if (value != 1) { + v.push_back(value % 2); + vhdl_int_to_unsigned(value / 2, v); + } else { + v.push_back(1); + } + return; +} + +void vhdl_int_to_signed(int value, std::vector &v) { + if (value < 0) { + int val = 1; + int size = 1; + while (val < -value) { + val *= 2; + size += 1; + } + vhdl_int_to_unsigned(val + value, v); + for (int i = v.size(); i < size - 1; i++) { + v.push_back(0); + } + v.push_back(1); + } else { + vhdl_int_to_unsigned(value, v); + v.push_back(0); + } + return; +} + +void vhdl_resize_unsigned(std::vector &v, int new_size) { + for (int i = v.size(); i < new_size; i++) { + v.push_back(0); + } +} + +void vhdl_resize_signed(std::vector &v, int new_size) { + int elem = 0; + if (v[v.size() - 1] == 1) + elem = 1; + for (int i = v.size(); i < new_size; i++) { + v.push_back(elem); + } +} + +bool vhdl_resize_signed_ok(std::vector v, int new_size) { + for (size_t i = v.size() - 1 - 1; i >= v.size() - 1 - (v.size() - new_size); i--) { + if (v[i] != v[v.size() - 1]) + return false; + } + return true; +}; + +bool vhdl_resize_unsigned_ok(std::vector v, int new_size) { + for (size_t i = v.size() - 1; i >= v.size() - 1 + 1 - (v.size() - new_size); i--) { + if (v[i] != 0) + return false; + } + return true; +};