diff --git a/DataFormats/L1DTTrackFinder/BuildFile.xml b/DataFormats/L1DTTrackFinder/BuildFile.xml index bd2a1d535a5e3..54b1c2a65f97d 100644 --- a/DataFormats/L1DTTrackFinder/BuildFile.xml +++ b/DataFormats/L1DTTrackFinder/BuildFile.xml @@ -3,4 +3,5 @@ + diff --git a/DataFormats/L1TMuonPhase2/interface/KMTFTrack.h b/DataFormats/L1TMuonPhase2/interface/KMTFTrack.h index 7e6197eea6d9e..a414b3458e169 100644 --- a/DataFormats/L1TMuonPhase2/interface/KMTFTrack.h +++ b/DataFormats/L1TMuonPhase2/interface/KMTFTrack.h @@ -10,30 +10,40 @@ namespace l1t { class KMTFTrack; typedef std::vector KMTFTrackCollection; typedef BXVector KMTFTrackBxCollection; + typedef math::Error<5>::type CovarianceMatrix5dim; + typedef math::Error<2>::type CovarianceMatrix2dim; class KMTFTrack : public reco::LeafCandidate { public: KMTFTrack() : reco::LeafCandidate(-1, reco::LeafCandidate::PolarLorentzVector(0.1, 0.0, 0.0, 0.105)), unconstrainedP4_(reco::LeafCandidate::PolarLorentzVector(0.1, 0.0, 0.0, 0.105)), - covariance_(std::vector(6, 0.0)), + covariance_(std::vector(15, 0.0)), + covarianceNB_(std::vector(3, 0.0)), curvVertex_(0), ptC_(0), phiVertex_(0), dxy_(0), + zVertex_(0), + kSlopeVertex_(0), curvMuon_(0), ptU_(0), phiMuon_(0), phiBMuon_(0), + zMuon_(0), + kSlopeMuon_(0), curv_(0), phi_(0), phiB_(0), + z_(0), + kSlope_(0), coarseEta_(0), approxPromptChi2_(0), approxPromptErrChi2_(0), approxDispChi2_(0), approxDispErrChi2_(0), hitPattern_(0), + thetaDigiPattern_(0), step_(1), sector_(0), wheel_(0), @@ -46,27 +56,35 @@ namespace l1t { ~KMTFTrack() override = default; - KMTFTrack(const l1t::MuonStubRef& seed, int phi, int phiB) + KMTFTrack(const l1t::MuonStubRef& seed, int phi, int phiB, int z, int kSlope) : reco::LeafCandidate(-1, reco::LeafCandidate::PolarLorentzVector(0.1, 0.0, 0.0, 0.105)), unconstrainedP4_(reco::LeafCandidate::PolarLorentzVector(0.1, 0.0, 0.0, 0.105)), - covariance_(std::vector(6, 0.0)), + covariance_(std::vector(15, 0.0)), + covarianceNB_(std::vector(3, 0.0)), curvVertex_(0), ptC_(0), phiVertex_(0), dxy_(0), + zVertex_(0), + kSlopeVertex_(0), curvMuon_(0), ptU_(0), phiMuon_(0), phiBMuon_(0), + zMuon_(0), + kSlopeMuon_(0), curv_(0), phi_(phi), phiB_(phiB), + z_(z), + kSlope_(kSlope), coarseEta_(0), approxPromptChi2_(0), approxPromptErrChi2_(0), approxDispChi2_(0), approxDispErrChi2_(0), hitPattern_(0), + thetaDigiPattern_(0), step_(seed->depthRegion()), sector_(seed->phiRegion()), wheel_(seed->etaRegion()), @@ -80,6 +98,7 @@ namespace l1t { residuals_.push_back(0); residuals_.push_back(0); residuals_.push_back(0); + residuals_.push_back(0); } reco::LeafCandidate::PolarLorentzVector displacedP4() const { return unconstrainedP4_; } @@ -92,6 +111,10 @@ namespace l1t { int phiAtMuon() const { return phiMuon_; } //unconstrained phiB at station 1 int phiBAtMuon() const { return phiBMuon_; } + //unconstrained z at station 1 + int zAtMuon() const { return zMuon_; } + //unconstrained kSLope at station 1 + int kSlopeAtMuon() const { return kSlopeMuon_; } //constrained pt int ptPrompt() const { return ptC_; } @@ -99,6 +122,10 @@ namespace l1t { int curvatureAtVertex() const { return curvVertex_; } //constrained phi at the vertex int phiAtVertex() const { return phiVertex_; } + //constrained z at the vertex + int zAtVertex() const { return zVertex_; } + //constrained kSlope at the vertex + int kSlopeAtVertex() const { return kSlopeVertex_; } //Impact parameter as calculated from the muon track int dxy() const { return dxy_; } //Unconstrained curvature at the Muon systen @@ -107,6 +134,10 @@ namespace l1t { int positionAngle() const { return phi_; } //Unconstrained bending angle at the Muon systen int bendingAngle() const { return phiB_; } + + //global z and slope of stub + int zPosition() const { return z_; } + int kSlope() const { return kSlope_; } //Coarse eta caluclated only using phi segments int coarseEta() const { return coarseEta_; } //Approximate Chi2 metrics @@ -116,6 +147,7 @@ namespace l1t { int approxDispErrChi2() const { return approxDispErrChi2_; } int hitPattern() const { return hitPattern_; } + int thetaDigiPattern() const { return thetaDigiPattern_; } //step; int step() const { return step_; } //sector; @@ -156,8 +188,37 @@ namespace l1t { return kalmanGain0_; } + const std::vector& ThetaGain1D(unsigned int step) const { + switch (step) { + case 3: + return ThetaGain1D3_; + case 2: + return ThetaGain1D2_; + case 1: + return ThetaGain1D1_; + case 0: + return ThetaGain1D0_; + } + return ThetaGain1D0_; + } + + const std::vector& ThetaGain(unsigned int step) const { + switch (step) { + case 3: + return ThetaGain3_; + case 2: + return ThetaGain2_; + case 1: + return ThetaGain1_; + case 0: + return ThetaGain0_; + } + return ThetaGain0_; + } + //get covariance const std::vector& covariance() const { return covariance_; } + const std::vector& covarianceNB() const { return covarianceNB_; } //get residual int residual(uint i) const { return residuals_[i]; } @@ -188,23 +249,29 @@ namespace l1t { } //Set coordinates general - void setCoordinates(int step, int curv, int phi, int phiB) { + void setCoordinates(int step, int curv, int phi, int phiB, int z, int kSlope) { step_ = step; curv_ = curv; phiB_ = phiB; phi_ = phi; + z_ = z; + kSlope_ = kSlope; } - void setCoordinatesAtVertex(int curv, int phi, int dxy) { + void setCoordinatesAtVertex(int curv, int phi, int dxy, int z, int kSlope) { curvVertex_ = curv; phiVertex_ = phi; dxy_ = dxy; + zVertex_ = z; + kSlopeVertex_ = kSlope; } - void setCoordinatesAtMuon(int curv, int phi, int phiB) { + void setCoordinatesAtMuon(int curv, int phi, int phiB, int z, int kSlope) { curvMuon_ = curv; phiMuon_ = phi; phiBMuon_ = phiB; + zMuon_ = z; + kSlopeMuon_ = kSlope; } void setPt(int ptC, int ptU) { @@ -215,6 +282,7 @@ namespace l1t { void setCoarseEta(int eta) { coarseEta_ = eta; } void setHitPattern(int pattern) { hitPattern_ = pattern; } + void setThetaDigiPattern(int theta_pattern) { thetaDigiPattern_ = theta_pattern; } void setApproxChi2(int chi, int chiErr, bool prompt) { if (prompt) { @@ -305,14 +373,101 @@ namespace l1t { } } + void setThetaGain1D(unsigned int step, + unsigned int K, + int priorThetaPattern, + int seedStation, + int priorPhiPattern, + float G31, + float G32, + float G41, + float G42) { + std::vector* v1 = nullptr; + switch (step) { + case 3: + v1 = &ThetaGain1D3_; + break; + case 2: + v1 = &ThetaGain1D2_; + break; + case 1: + v1 = &ThetaGain1D1_; + break; + case 0: + v1 = &ThetaGain1D0_; + break; + default: + throw cms::Exception("WrongCondition") << "1D: Critical ERROR on setting the Theta gain\n"; + } + v1->push_back(static_cast(K)); + v1->push_back(static_cast(priorThetaPattern)); + v1->push_back(static_cast(seedStation)); + v1->push_back(static_cast(priorPhiPattern)); + v1->push_back(G31); + v1->push_back(G32); + v1->push_back(G41); + v1->push_back(G42); + } + + void setThetaGain(unsigned int step, + unsigned int K, + int priorThetaPattern, + int seedStation, + int priorPhiPattern, + float G32, + float G33, + float G42, + float G43) { + std::vector* v2 = nullptr; + switch (step) { + case 3: + v2 = &ThetaGain3_; + break; + case 2: + v2 = &ThetaGain2_; + break; + case 1: + v2 = &ThetaGain1_; + break; + case 0: + v2 = &ThetaGain0_; + break; + default: + throw cms::Exception("WrongCondition") << "Critical ERROR on setting the Theta gain\n"; + } + v2->push_back(static_cast(K)); + v2->push_back(static_cast(priorThetaPattern)); + v2->push_back(static_cast(seedStation)); + v2->push_back(static_cast(priorPhiPattern)); + v2->push_back(G32); + v2->push_back(G33); + v2->push_back(G42); + v2->push_back(G43); + } + //set covariance - void setCovariance(const CovarianceMatrix& c) { + void setCovariance(const CovarianceMatrix5dim& c) { covariance_[0] = c(0, 0); covariance_[1] = c(0, 1); covariance_[2] = c(1, 1); covariance_[3] = c(0, 2); covariance_[4] = c(1, 2); covariance_[5] = c(2, 2); + covariance_[6] = c(0, 3); + covariance_[7] = c(1, 3); + covariance_[8] = c(2, 3); + covariance_[9] = c(3, 3); + covariance_[10] = c(0, 4); + covariance_[11] = c(1, 4); + covariance_[12] = c(2, 4); + covariance_[13] = c(3, 4); + covariance_[14] = c(4, 4); + } + + void setCovarianceNB(const CovarianceMatrix2dim& c) { + covarianceNB_[0] = c(0, 0); + covarianceNB_[1] = c(0, 1); + covarianceNB_[2] = c(1, 1); } //set fine eta @@ -329,6 +484,7 @@ namespace l1t { //Covariance matrix for studies std::vector covariance_; + std::vector covarianceNB_; l1t::MuonStubRefVector stubs_; //vertex coordinates @@ -336,17 +492,23 @@ namespace l1t { int ptC_; int phiVertex_; int dxy_; + int zVertex_; + int kSlopeVertex_; //muon coordinates int curvMuon_; int ptU_; int phiMuon_; int phiBMuon_; + int zMuon_; + int kSlopeMuon_; //generic coordinates int curv_; int phi_; int phiB_; + int z_; + int kSlope_; //common coordinates int coarseEta_; @@ -358,6 +520,8 @@ namespace l1t { //phi bitmask int hitPattern_; + //bitmask pattern based on theta digi presence + int thetaDigiPattern_; //propagation step int step_; @@ -393,6 +557,16 @@ namespace l1t { std::vector kalmanGain3_; std::vector residuals_; + + std::vector ThetaGain1D0_; + std::vector ThetaGain1D1_; + std::vector ThetaGain1D2_; + std::vector ThetaGain1D3_; + + std::vector ThetaGain0_; + std::vector ThetaGain1_; + std::vector ThetaGain2_; + std::vector ThetaGain3_; }; } // namespace l1t diff --git a/DataFormats/L1TMuonPhase2/src/classes_def.xml b/DataFormats/L1TMuonPhase2/src/classes_def.xml index 4c34f5483adcf..45b0a62ee6246 100644 --- a/DataFormats/L1TMuonPhase2/src/classes_def.xml +++ b/DataFormats/L1TMuonPhase2/src/classes_def.xml @@ -22,7 +22,8 @@ - + + diff --git a/L1Trigger/Phase2L1GMT/interface/KMTFCore.h b/L1Trigger/Phase2L1GMT/interface/KMTFCore.h index 19291d1f53031..99cb6553f3522 100644 --- a/L1Trigger/Phase2L1GMT/interface/KMTFCore.h +++ b/L1Trigger/Phase2L1GMT/interface/KMTFCore.h @@ -2,6 +2,8 @@ Kalman Filter L1 Muon algorithm Tyler Lam (UCLA) Sep. 2021 +Delano Campos (UCLA) +Feb. 2026 */ #ifndef L1Trigger_Phase2L1GMT_KMTFCore_h #define L1Trigger_Phase2L1GMT_KMTFCore_h @@ -16,12 +18,25 @@ namespace Phase2L1GMT { class KMTFCore { public: typedef ROOT::Math::SVector Vector2; + typedef ROOT::Math::SVector Vector3; + typedef ROOT::Math::SVector Vector4; typedef ROOT::Math::SMatrix > CovarianceMatrix2; + typedef ROOT::Math::SMatrix > CovarianceMatrix3; + typedef ROOT::Math::SMatrix > CovarianceMatrix4; typedef ROOT::Math::SMatrix Matrix32; + typedef ROOT::Math::SMatrix Matrix53; typedef ROOT::Math::SMatrix Matrix23; + typedef ROOT::Math::SMatrix Matrix45; + typedef ROOT::Math::SMatrix Matrix54; typedef ROOT::Math::SMatrix Matrix13; + typedef ROOT::Math::SMatrix Matrix35; typedef ROOT::Math::SMatrix Matrix31; typedef ROOT::Math::SMatrix Matrix33; + typedef ROOT::Math::SMatrix Matrix55; + typedef ROOT::Math::SMatrix Matrix15; + typedef ROOT::Math::SMatrix Matrix51; + typedef ROOT::Math::SMatrix Matrix52; + typedef ROOT::Math::SMatrix Matrix25; KMTFCore(const edm::ParameterSet& settings); @@ -44,12 +59,14 @@ namespace Phase2L1GMT { void vertexConstraintOffline(l1t::KMTFTrack& track); void vertexConstraintLUT(l1t::KMTFTrack& track); int hitPattern(const l1t::KMTFTrack& track); + int thetaDigiPattern(const l1t::KMTFTrack& track); int customBitmask(unsigned int bit1, unsigned int bit2, unsigned int bit3, unsigned int bit4); bool getBit(int bitmask, int pos); void setFourVectors(l1t::KMTFTrack& track); bool estimateChiSquare(l1t::KMTFTrack& track, bool vertex); void setRank(l1t::KMTFTrack& track, bool vertex); int wrapAround(int value, int maximum); + int satKSlope(int k); int encode(bool ownwheel, int sector, int tag); std::pair getByCode(const std::vector& tracks, int mask); uint twosCompToBits(int q); @@ -115,9 +132,13 @@ namespace Phase2L1GMT { //bits for fixed point precision static const int PHIBSCALE = 16; static const int PHIBSCALE_INT = 5; + static const int ZDELTAR_BITS = 14; + static const int ZDELTAR_BITSINT = 1; static const int BITSCURV = 16; static const int BITSPHI = 18; static const int BITSPHIB = 17; // 12 bits *28 (+5 bits) + static const int BITSZ = 16; + static const int BITSKSLOPE = 16; static const int BITSPARAM = 14; static const int GAIN_0 = 9; static const int GAIN_0INT = 6; @@ -134,12 +155,26 @@ namespace Phase2L1GMT { static const int GAIN2_4INT = 4; static const int GAIN2_5 = 12; static const int GAIN2_5INT = 0; + + static const int GAINT_6 = 12; //new theta view gain bits. + static const int GAINT_6INT = 3; + static const int GAINT_7 = 12; + static const int GAINT_7INT = 3; + static const int GAINT_8 = 12; + static const int GAINT_8INT = 3; + static const int GAINT_9 = 12; + static const int GAINT_9INT = 3; //STUFF NOT USED IN THE FIRMWARE BUT ONLY FOR DEBUGGING /////////////////////////////////////////////////////// bool useOfflineAlgo_; + //used for z propagation + std::vector zdeltaR_; + //multiple scattering constants used in multiple scattering matrix to build RMS values std::vector mScatteringPhi_; std::vector mScatteringPhiB_; + std::vector mScatteringz_; + std::vector mScatteringkSlope_; //point resolution for phi double pointResolutionPhi_; //point resolution for phiB @@ -149,6 +184,10 @@ namespace Phase2L1GMT { //double pointResolutionPhiB_; //point resolution for vertex double pointResolutionVertex_; + // point resolution for z + std::vector pointResolutionz_; + // point resolution for kSlope + std::vector pointResolutionkSlope_; std::vector curvResolution1_; std::vector curvResolution2_; //Sorter diff --git a/L1Trigger/Phase2L1GMT/interface/KMTFLUTs.h b/L1Trigger/Phase2L1GMT/interface/KMTFLUTs.h index 42027a4c38e43..28a1f99c233a3 100644 --- a/L1Trigger/Phase2L1GMT/interface/KMTFLUTs.h +++ b/L1Trigger/Phase2L1GMT/interface/KMTFLUTs.h @@ -6,69 +6,94 @@ #include "TFile.h" #include #include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/Utilities/interface/Exception.h" namespace Phase2L1GMT { class KMTFLUTs { public: - KMTFLUTs(const std::string &filename) { + KMTFLUTs(const std::string &filename, const std::string &ThetaFilename) { edm::FileInPath path(filename); - lutFile_ = new TFile(path.fullPath().c_str()); - lut_[3 * 64 + 8] = (TH1 *)lutFile_->Get("gain_8_3"); - lut_[2 * 64 + 8] = (TH1 *)lutFile_->Get("gain_8_2"); - lut_[2 * 64 + 12] = (TH1 *)lutFile_->Get("gain_12_2"); - lut_[2 * 64 + 4] = (TH1 *)lutFile_->Get("gain_4_2"); - lut_[1 * 64 + 12] = (TH1 *)lutFile_->Get("gain_12_1"); - lut_[1 * 64 + 10] = (TH1 *)lutFile_->Get("gain_10_1"); - lut_[1 * 64 + 6] = (TH1 *)lutFile_->Get("gain_6_1"); - lut_[1 * 64 + 14] = (TH1 *)lutFile_->Get("gain_14_1"); - lut_[3] = (TH1 *)lutFile_->Get("gain_3_0"); - lut_[5] = (TH1 *)lutFile_->Get("gain_5_0"); - lut_[6] = (TH1 *)lutFile_->Get("gain_6_0"); - lut_[7] = (TH1 *)lutFile_->Get("gain_7_0"); - lut_[9] = (TH1 *)lutFile_->Get("gain_9_0"); - lut_[10] = (TH1 *)lutFile_->Get("gain_10_0"); - lut_[11] = (TH1 *)lutFile_->Get("gain_11_0"); - lut_[12] = (TH1 *)lutFile_->Get("gain_12_0"); - lut_[13] = (TH1 *)lutFile_->Get("gain_13_0"); - lut_[14] = (TH1 *)lutFile_->Get("gain_14_0"); - lut_[15] = (TH1 *)lutFile_->Get("gain_15_0"); - - lut2HH_[3 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_3_HH"); - lut2HH_[2 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_2_HH"); - lut2HH_[2 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_2_HH"); - lut2HH_[1 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_1_HH"); - lut2HH_[1 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_1_HH"); - lut2HH_[1 * 64 + 2] = (TH1 *)lutFile_->Get("gain2_2_1_HH"); - - lut2LH_[3 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_3_LH"); - lut2LH_[2 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_2_LH"); - lut2LH_[2 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_2_LH"); - lut2LH_[1 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_1_LH"); - lut2LH_[1 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_1_LH"); - lut2LH_[1 * 64 + 2] = (TH1 *)lutFile_->Get("gain2_2_1_LH"); - - lut2HL_[3 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_3_HL"); - lut2HL_[2 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_2_HL"); - lut2HL_[2 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_2_HL"); - lut2HL_[1 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_1_HL"); - lut2HL_[1 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_1_HL"); - lut2HL_[1 * 64 + 2] = (TH1 *)lutFile_->Get("gain2_2_1_HL"); - - lut2LL_[3 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_3_LL"); - lut2LL_[2 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_2_LL"); - lut2LL_[2 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_2_LL"); - lut2LL_[1 * 64 + 8] = (TH1 *)lutFile_->Get("gain2_8_1_LL"); - lut2LL_[1 * 64 + 4] = (TH1 *)lutFile_->Get("gain2_4_1_LL"); - lut2LL_[1 * 64 + 2] = (TH1 *)lutFile_->Get("gain2_2_1_LL"); - - coarseEta_ = (TH1 *)lutFile_->Get("coarseETALUT"); + edm::FileInPath pathTheta(ThetaFilename); + lutFile_ = std::unique_ptr{TFile::Open(path.fullPath().c_str())}; + if (!lutFile_ || lutFile_->IsZombie()) { + throw cms::Exception("KMTFLUT") << "Missing packedGainLUTs file\n"; + } + lutThetaFile_ = std::unique_ptr(TFile::Open(pathTheta.fullPath().c_str())); + if (!lutThetaFile_ || lutThetaFile_->IsZombie()) { + throw cms::Exception("KMTFLUT") << "Missing packedThetaGainLUTs file\n"; + }; + lut_[3 * 64 + 8] = static_cast(lutFile_->Get("gain_8_3")); + lut_[2 * 64 + 8] = static_cast(lutFile_->Get("gain_8_2")); + lut_[2 * 64 + 12] = static_cast(lutFile_->Get("gain_12_2")); + lut_[2 * 64 + 4] = static_cast(lutFile_->Get("gain_4_2")); + lut_[1 * 64 + 12] = static_cast(lutFile_->Get("gain_12_1")); + lut_[1 * 64 + 10] = static_cast(lutFile_->Get("gain_10_1")); + lut_[1 * 64 + 6] = static_cast(lutFile_->Get("gain_6_1")); + lut_[1 * 64 + 14] = static_cast(lutFile_->Get("gain_14_1")); + lut_[3] = static_cast(lutFile_->Get("gain_3_0")); + lut_[5] = static_cast(lutFile_->Get("gain_5_0")); + lut_[6] = static_cast(lutFile_->Get("gain_6_0")); + lut_[7] = static_cast(lutFile_->Get("gain_7_0")); + lut_[9] = static_cast(lutFile_->Get("gain_9_0")); + lut_[10] = static_cast(lutFile_->Get("gain_10_0")); + lut_[11] = static_cast(lutFile_->Get("gain_11_0")); + lut_[12] = static_cast(lutFile_->Get("gain_12_0")); + lut_[13] = static_cast(lutFile_->Get("gain_13_0")); + lut_[14] = static_cast(lutFile_->Get("gain_14_0")); + lut_[15] = static_cast(lutFile_->Get("gain_15_0")); + + lut2HH_[3 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_3_HH")); + lut2HH_[2 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_2_HH")); + lut2HH_[2 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_2_HH")); + lut2HH_[1 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_1_HH")); + lut2HH_[1 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_1_HH")); + lut2HH_[1 * 64 + 2] = static_cast(lutFile_->Get("gain2_2_1_HH")); + + lut2LH_[3 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_3_LH")); + lut2LH_[2 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_2_LH")); + lut2LH_[2 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_2_LH")); + lut2LH_[1 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_1_LH")); + lut2LH_[1 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_1_LH")); + lut2LH_[1 * 64 + 2] = static_cast(lutFile_->Get("gain2_2_1_LH")); + + lut2HL_[3 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_3_HL")); + lut2HL_[2 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_2_HL")); + lut2HL_[2 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_2_HL")); + lut2HL_[1 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_1_HL")); + lut2HL_[1 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_1_HL")); + lut2HL_[1 * 64 + 2] = static_cast(lutFile_->Get("gain2_2_1_HL")); + + lut2LL_[3 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_3_LL")); + lut2LL_[2 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_2_LL")); + lut2LL_[2 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_2_LL")); + lut2LL_[1 * 64 + 8] = static_cast(lutFile_->Get("gain2_8_1_LL")); + lut2LL_[1 * 64 + 4] = static_cast(lutFile_->Get("gain2_4_1_LL")); + lut2LL_[1 * 64 + 2] = static_cast(lutFile_->Get("gain2_2_1_LL")); + + coarseEta_ = static_cast(lutFile_->Get("coarseETALUT")); + + lutTheta1D_[1 * 64 + 0] = static_cast(lutThetaFile_->Get("gain1D_0_1")); + lutTheta1D_[1 * 64 + 2] = static_cast(lutThetaFile_->Get("gain1D_2_1")); + lutTheta1D_[1 * 64 + 4] = static_cast(lutThetaFile_->Get("gain1D_4_1")); + lutTheta1D_[1 * 64 + 6] = static_cast(lutThetaFile_->Get("gain1D_6_1")); + lutTheta1D_[2 * 64 + 4] = static_cast(lutThetaFile_->Get("gain1D_4_2")); + lutTheta1D_[3 * 64 + 0] = static_cast(lutThetaFile_->Get("gain1D_0_3")); + lutTheta2D_[1 * 64 + 0] = static_cast(lutThetaFile_->Get("gain2D_0_1")); + lutTheta2D_[1 * 64 + 2] = static_cast(lutThetaFile_->Get("gain2D_2_1")); + lutTheta2D_[1 * 64 + 4] = static_cast(lutThetaFile_->Get("gain2D_4_1")); + lutTheta2D_[1 * 64 + 6] = static_cast(lutThetaFile_->Get("gain2D_6_1")); + lutTheta2D_[2 * 64 + 0] = static_cast(lutThetaFile_->Get("gain2D_0_2")); + lutTheta2D_[2 * 64 + 4] = static_cast(lutThetaFile_->Get("gain2D_4_2")); + lutTheta2D_[3 * 64 + 0] = static_cast(lutThetaFile_->Get("gain2D_0_3")); + lutTheta1D11_[2 * 64 + 0] = static_cast(lutThetaFile_->Get("gain1D_0_2_phi1100")); + lutTheta1D01_[2 * 64 + 0] = static_cast(lutThetaFile_->Get("gain1D_0_2_phi0100")); + lutTheta1D10_[2 * 64 + 0] = static_cast(lutThetaFile_->Get("gain1D_0_2_phi1000")); } ~KMTFLUTs() { lutFile_->Close(); - if (lutFile_ != nullptr) - delete lutFile_; + lutThetaFile_->Close(); } std::vector trackGain(uint step, uint bitmask, uint K) { @@ -110,12 +135,51 @@ namespace Phase2L1GMT { return uint((1 << 12) * coarseEta_->GetBinContent(coarseEta_->GetXaxis()->FindBin(mask)) / M_PI); } - TFile *lutFile_; + std::vector trackGainTheta(uint step, uint bitmask, uint K, bool is2D) const { + std::vector gain(4, 0.0); + const TH1 *h; + if (is2D) { + h = lutTheta2D_.at(64 * step + bitmask); + } else { + h = lutTheta1D_.at(64 * step + bitmask); + } + gain[0] = h->GetBinContent(K + 1); + gain[1] = h->GetBinContent(512 + K + 1); + gain[2] = h->GetBinContent(2 * 512 + K + 1); + gain[3] = h->GetBinContent(3 * 512 + K + 1); + return gain; + } + + std::vector trackGainTheta2(uint step, uint bitmask, uint phiBitmask, uint K) const { + std::vector gain(4, 0.0); + const TH1 *h; + if (phiBitmask == 0b1100) { + h = lutTheta1D11_.at(64 * step + bitmask); + } else if (phiBitmask == 0b1000) { + h = lutTheta1D10_.at(64 * step + bitmask); + } else { + h = lutTheta1D01_.at(64 * step + bitmask); + } + gain[0] = h->GetBinContent(K + 1); + gain[1] = h->GetBinContent(512 + K + 1); + gain[2] = h->GetBinContent(2 * 512 + K + 1); + gain[3] = h->GetBinContent(3 * 512 + K + 1); + return gain; + } + + std::unique_ptr lutFile_; + std::unique_ptr lutThetaFile_; std::map lut_; std::map lut2HH_; std::map lut2LH_; std::map lut2HL_; std::map lut2LL_; + + std::map lutTheta1D_; + std::map lutTheta2D_; + std::map lutTheta1D11_; + std::map lutTheta1D10_; + std::map lutTheta1D01_; const TH1 *coarseEta_; }; diff --git a/L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTBarrelStubProcessor.h b/L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTBarrelStubProcessor.h index 81607c571c702..523ef4649e865 100644 --- a/L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTBarrelStubProcessor.h +++ b/L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTBarrelStubProcessor.h @@ -5,6 +5,8 @@ #include "DataFormats/L1DTTrackFinder/interface/L1Phase2MuDTPhContainer.h" #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambThDigi.h" #include "DataFormats/L1DTTrackFinder/interface/L1MuDTChambThContainer.h" +#include "DataFormats/L1DTTrackFinder/interface/L1Phase2MuDTExtPhiThetaPair.h" +#include "DataFormats/L1DTTrackFinder/interface/L1Phase2MuDTExtPhiThetaPairContainer.h" #include "DataFormats/L1TMuonPhase2/interface/MuonStub.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "CondFormats/L1TObjects/interface/L1TMuonBarrelParams.h" @@ -17,9 +19,10 @@ class L1TPhase2GMTBarrelStubProcessor { ~L1TPhase2GMTBarrelStubProcessor(); - l1t::MuonStubCollection makeStubs(const L1Phase2MuDTPhContainer*, const L1MuDTChambThContainer*); + l1t::MuonStubCollection makeStubs(const L1Phase2MuDTExtPhiThetaPairContainer*); private: + l1t::MuonStub buildStubwithZandkSlope(const L1Phase2MuDTExtPhiThetaPair&); l1t::MuonStub buildStub(const L1Phase2MuDTPhDigi&, const L1MuDTChambThDigi*); l1t::MuonStub buildStubNoEta(const L1Phase2MuDTPhDigi&); diff --git a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTKMTFProducer.cc b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTKMTFProducer.cc index 04855cd154cb3..c8c3a6b0411fd 100644 --- a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTKMTFProducer.cc +++ b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTKMTFProducer.cc @@ -95,12 +95,18 @@ void Phase2L1TGMTKMTFProducer::produce(edm::Event& iEvent, const edm::EventSetup std::vector kmtfTracks; for (const auto& track : kmtfOutput.first) { kmtfTracks.push_back(track); + //convert kmtf track.zPosition (LSB_kmtf=1500./2<<16) into z0 for SAMuon. + //LSBSAz0 = 1.6cm is too small given only 5 bits. the range fro z0 is +- 25.6cm from DataFormats/L1TMuonPhase2/interface/Constants.h + //eff. LSB = 1.6cm * 16 = 25.6 cm to widen the z0 range that can covered given only 5 bits. + //full conversion is z0_SA=z0_kmtf*(LSB_kmtf)*(1/LSB_SA)*(1/16) to convert between kmtf->SAMuon digitizations + ap_fixed z0 = + track.zPosition() * ap_ufixed<18, 0>((1500.0 / (1 << 16)) / (LSBSAz0 * 16)); l1t::SAMuon p(track.p4(), (track.curvatureAtVertex() < 0), track.ptPrompt(), - track.coarseEta(), + round(track.eta() / (M_PI / (1 << 12))), track.phiAtMuon() / (1 << 5), - 0, + z0, 0, track.stubs().size() - 1); p.setTF(l1t::tftype::bmtf); @@ -123,12 +129,14 @@ void Phase2L1TGMTKMTFProducer::produce(edm::Event& iEvent, const edm::EventSetup for (const auto& track : kmtfOutput.second) { kmtfTracks.push_back(track); ap_int<7> dxy = track.dxy() * ap_ufixed<8, 1>(1.606); + ap_fixed z0 = + track.zPosition() * ap_ufixed<18, 0>((1500.0 / (1 << 16)) / (LSBSAz0 * 16)); l1t::SAMuon p(track.displacedP4(), (track.curvatureAtMuon() < 0), track.ptDisplaced(), - track.coarseEta(), + round(track.eta() / (M_PI / (1 << 12))), track.phiAtMuon() / (1 << 5), - 0, + z0, dxy, track.approxDispChi2()); p.setTF(l1t::tftype::bmtf); diff --git a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTStubProducer.cc b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTStubProducer.cc index ed680955534ea..3f65c2aacd3db 100644 --- a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTStubProducer.cc +++ b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTStubProducer.cc @@ -11,6 +11,8 @@ #include "L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTEndcapStubProcessor.h" #include "L1Trigger/Phase2L1GMT/interface/L1TPhase2GMTBarrelStubProcessor.h" +#include "DataFormats/L1DTTrackFinder/interface/L1Phase2MuDTExtPhiThetaPairContainer.h" + #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/ConsumesCollector.h" #include "FWCore/Framework/interface/ESProducts.h" @@ -34,6 +36,7 @@ class Phase2L1TGMTStubProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT> srcCSC_; edm::EDGetTokenT srcDT_; edm::EDGetTokenT srcDTTheta_; + edm::EDGetTokenT srcDTPairs_; edm::EDGetTokenT srcRPC_; L1TPhase2GMTEndcapStubProcessor* procEndcap_; @@ -47,6 +50,7 @@ Phase2L1TGMTStubProducer::Phase2L1TGMTStubProducer(const edm::ParameterSet& iCon consumes>(iConfig.getParameter("srcCSC"))), srcDT_(consumes(iConfig.getParameter("srcDT"))), srcDTTheta_(consumes(iConfig.getParameter("srcDTTheta"))), + srcDTPairs_(consumes(iConfig.getParameter("srcDTPairs"))), srcRPC_(consumes(iConfig.getParameter("srcRPC"))), procEndcap_(new L1TPhase2GMTEndcapStubProcessor(iConfig.getParameter("Endcap"))), procBarrel_(new L1TPhase2GMTBarrelStubProcessor(iConfig.getParameter("Barrel"))), @@ -107,6 +111,9 @@ void Phase2L1TGMTStubProducer::produce(edm::Event& iEvent, const edm::EventSetup Handle dtThetaDigis; iEvent.getByToken(srcDTTheta_, dtThetaDigis); + Handle dtPairs; + iEvent.getByToken(srcDTPairs_, dtPairs); + //Generate a unique stub ID l1t::MuonStubCollection stubs; l1t::MuonStubCollection stubsKMTF; @@ -115,7 +122,7 @@ void Phase2L1TGMTStubProducer::produce(edm::Event& iEvent, const edm::EventSetup for (auto& stub : stubsEndcap) { stubs.push_back(stub); } - l1t::MuonStubCollection stubsBarrel = procBarrel_->makeStubs(dtDigis.product(), dtThetaDigis.product()); + l1t::MuonStubCollection stubsBarrel = procBarrel_->makeStubs(dtPairs.product()); for (auto& stub : stubsBarrel) { //convert to Hybrid stubs.push_back(convertToHybrid(stub)); @@ -139,6 +146,7 @@ void Phase2L1TGMTStubProducer::fillDescriptions(edm::ConfigurationDescriptions& desc.add("srcCSC", edm::InputTag("simCscTriggerPrimitiveDigis")); desc.add("srcDT", edm::InputTag("dtTriggerPhase2PrimitiveDigis")); desc.add("srcDTTheta", edm::InputTag("simDtTriggerPrimitiveDigis")); + desc.add("srcDTPairs", edm::InputTag("dtTriggerPhase2PrimitivePairDigis")); desc.add("srcRPC", edm::InputTag("simMuonRPCDigis")); { edm::ParameterSetDescription psd0; diff --git a/L1Trigger/Phase2L1GMT/python/gmtKMTFMuons_cfi.py b/L1Trigger/Phase2L1GMT/python/gmtKMTFMuons_cfi.py index b328e18adac84..5119381b7c8fa 100644 --- a/L1Trigger/Phase2L1GMT/python/gmtKMTFMuons_cfi.py +++ b/L1Trigger/Phase2L1GMT/python/gmtKMTFMuons_cfi.py @@ -7,12 +7,14 @@ algo = cms.PSet( verbose = cms.bool(False), lutFile = cms.string("L1Trigger/Phase2L1GMT/data/packedGainLUTs.root"), + lutThetaFile = cms.string("L1Trigger/Phase2L1GMT/data/packedThetaGainLUTs.root"), initialK = cms.vdouble(-0.4576229536749278, -0.6364802777566145, -1.0305030909883524, -1.7272067322624118), initialK2 = cms.vdouble(-6.442002637356136e-05, -9.582709649965545e-05, -0.0002741064246218815, -0.0014910074450869175), eLoss = cms.vdouble(6.77765e-05,0,0,0), aPhi = cms.vdouble(5.6533349391874275, 0.03477876333443834, 0.032506522838098864, 0.024752809174909053), aPhiB = cms.vdouble(-2.02, -0.2994087741381382, -0.4033668521165302, -0.3592231728688621), + zdeltaR = cms.vdouble(445.0, 526.0-445.0, 635.0-526.0, 730.0-635.0), aPhiBNLO = cms.vdouble(9.04133e-05,0,0,0), bPhi = cms.vdouble(-1,.18245,.20898,.17286), @@ -58,11 +60,15 @@ ###Only for the offline algo -not in firmware -------------------- mScatteringPhi = cms.vdouble(0.1169021113298598, 0.00016777763395543814, 0.0004322078772344548, 0.00024985881710722107), mScatteringPhiB = cms.vdouble(.0522762, 0.01762000062188365, 0.03508319015441297, 0.03126825551530328), + mScatteringz = cms.vdouble(0.030625, 0.00010609, 0.00033489, 0.00033489), + mScatteringkSlope = cms.vdouble(0.083521, 0.032041, 0.047961, 0.047961), pointResolutionPhi = cms.double(1.), pointResolutionPhiB = cms.double(12493.7429036), pointResolutionPhiBH = cms.vdouble(19925.62953113343, 15583.06791339368, 10258.11768352221, 15462.112839170433), pointResolutionPhiBL = cms.vdouble(161519.85395846734, 155051.58394241595, 149693.88179343447, 174896.46766622085), pointResolutionVertex = cms.double(1.), + pointResolutionz = cms.vdouble(214369, 9980.01, 9254.44, 9254.44), + pointResolutionkSlope = cms.vdouble(186624, 350464, 375769, 375769), curvResolution1 = cms.vdouble(1, 2.36097e+03, 8.73003e+02, 2.58138e5), curvResolution2 = cms.vdouble(1, 4.903692e+00, 4.87941e+01, 0)), Nprompt = cms.uint32(12), diff --git a/L1Trigger/Phase2L1GMT/python/gmtStubs_cfi.py b/L1Trigger/Phase2L1GMT/python/gmtStubs_cfi.py index 1e4f32896bc7f..97053d7bcc162 100644 --- a/L1Trigger/Phase2L1GMT/python/gmtStubs_cfi.py +++ b/L1Trigger/Phase2L1GMT/python/gmtStubs_cfi.py @@ -5,6 +5,7 @@ srcCSC = cms.InputTag("simCscTriggerPrimitiveDigis"), srcDT = cms.InputTag("dtTriggerPhase2PrimitiveDigis"), srcDTTheta = cms.InputTag("simDtTriggerPrimitiveDigis"), + srcDTPairs = cms.InputTag("dtTriggerPhase2PrimitivePairDigis"), srcRPC = cms.InputTag("simMuonRPCDigis"), Endcap =cms.PSet( verbose = cms.uint32(0), diff --git a/L1Trigger/Phase2L1GMT/src/KMTF.cc b/L1Trigger/Phase2L1GMT/src/KMTF.cc index 918554895b131..b96989ee9cecb 100644 --- a/L1Trigger/Phase2L1GMT/src/KMTF.cc +++ b/L1Trigger/Phase2L1GMT/src/KMTF.cc @@ -72,74 +72,80 @@ std::pair, std::vector > KMTF::proce } bool pre_patterns = (verbose_ > 1) && ((Nstubs4 + Nstubs3 + Nstubs2 + Nstubs1) > 2); - //OK now process the data almost as in hardware for (unsigned int i = 0; i < 32; ++i) { //print the stubs taking into account - bool patterns = pre_patterns && ((i < Nstubs4) || (i < Nstubs3) || (i << Nstubs2)); + bool patterns = pre_patterns && ((i < Nstubs4) || (i < Nstubs3) || (i < Nstubs2)); if (patterns) { edm::LogInfo("KMTF") << "KMTFPattern " << std::flush; if (i < Nstubs4) edm::LogInfo("KMTF") << stubs4[0]->coord1() << " " << stubs4[0]->coord2() << " " << stubs4[0]->quality() << " " - << " 1 " << stubs4[0]->kmtf_address() << " 0 " << std::flush; + << " 1 " << stubs4[0]->kmtf_address() << " 0 " << stubs4[0]->eta1() << " " + << stubs4[0]->eta2() << " " << (stubs4[0]->etaQuality() > 0 ? 1 : 0) << " " << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; if (i < Nstubs3) { for (const auto& s : stubs3) { edm::LogInfo("KMTF") << "" << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " - << s->kmtf_address() << " 0 " << std::flush; + << s->kmtf_address() << " 0 " << s->eta1() << " " << s->eta2() << " " + << (s->etaQuality() > 0 ? 1 : 0) << " " << std::flush; } //pad with zeros for (unsigned int j = stubs3.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } } else { for (unsigned int j = stubs3.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } for (const auto& s : stubs3) { edm::LogInfo("KMTF") << "" << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " - << s->kmtf_address() << " 0 " << std::flush; + << s->kmtf_address() << " 0 " << s->eta1() << " " << s->eta2() << " " + << (s->etaQuality() > 0 ? 1 : 0) << " " << std::flush; } } if (i < Nstubs2) { for (const auto& s : stubs2) { edm::LogInfo("KMTF") << "" << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " - << s->kmtf_address() << " 0 " << std::flush; + << s->kmtf_address() << " 0 " << s->eta1() << " " << s->eta2() << " " + << (s->etaQuality() > 0 ? 1 : 0) << " " << std::flush; } //pad with zeros for (unsigned int j = stubs2.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } } else { for (unsigned int j = stubs2.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } for (const auto& s : stubs2) { edm::LogInfo("KMTF") << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " << s->kmtf_address() - << " 0 " << std::flush; + << " 0 " << s->eta1() << " " << s->eta2() << " " << (s->etaQuality() > 0 ? 1 : 0) << " " + << std::flush; } } if (i < Nstubs1) { for (const auto& s : stubs1) { edm::LogInfo("KMTF") << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " << s->kmtf_address() - << " 0 " << std::flush; + << " 0 " << s->eta1() << " " << s->eta2() << " " << (s->etaQuality() > 0 ? 1 : 0) << " " + << std::flush; } //pad with zeros for (unsigned int j = stubs1.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } } else { for (unsigned int j = stubs1.size(); j < 32; ++j) { - edm::LogInfo("KMTF") << "0 0 0 0 511 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 511 0 0 0 0 " << std::flush; } for (const auto& s : stubs1) { edm::LogInfo("KMTF") << s->coord1() << " " << s->coord2() << " " << s->quality() << " 1 " << s->kmtf_address() - << " 0 " << std::flush; + << " 0 " << s->eta1() << " " << s->eta2() << " " << (s->etaQuality() > 0 ? 1 : 0) << " " + << std::flush; } } } @@ -162,23 +168,24 @@ std::pair, std::vector > KMTF::proce if (tracks.first.id() & 0x1) edm::LogInfo("KMTF") << "1 " << (tracks.first.curvatureAtVertex() < 0 ? 1 : 0) << " " << tracks.first.ptPrompt() << " " << tracks.first.phiAtMuon() / (1 << 5) << " " - << tracks.first.coarseEta() << " " << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.first.rankPrompt() << " " << std::flush; - + << int(round(tracks.first.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.first.rankPrompt() + << " " << tracks.first.zPosition() << " " << tracks.first.kSlope() << " " << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; if (tracks.second.id() & 0x2) - edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtVertex() < 0 ? 1 : 0) << " " + edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtMuon() < 0 ? 1 : 0) << " " << tracks.second.ptDisplaced() << " " << tracks.second.phiAtMuon() / (1 << 5) << " " - << tracks.second.coarseEta() << " " << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.second.rankDisp() << " " << std::flush; - + << int(round(tracks.second.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.second.rankDisp() + << " " << tracks.second.zPosition() << " " << tracks.second.kSlope() << " " + << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; } } else if (patterns) { - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; } if (i < Nstubs3) { @@ -196,23 +203,26 @@ std::pair, std::vector > KMTF::proce if (tracks.first.id() & 0x1) edm::LogInfo("KMTF") << "1 " << (tracks.first.curvatureAtVertex() < 0 ? 1 : 0) << " " << tracks.first.ptPrompt() << " " << tracks.first.phiAtMuon() / (1 << 5) << " " - << tracks.first.coarseEta() << " " << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.first.rankPrompt() << " " << std::flush; + << int(round(tracks.first.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.first.rankPrompt() + << " " << tracks.first.zPosition() << " " << tracks.first.kSlope() << " " << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; if (tracks.second.id() & 0x2) - edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtVertex() < 0 ? 1 : 0) << " " + edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtMuon() < 0 ? 1 : 0) << " " << tracks.second.ptDisplaced() << " " << tracks.second.phiAtMuon() / (1 << 5) << " " - << tracks.second.coarseEta() << " " << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.second.rankDisp() << " " << std::flush; + << int(round(tracks.second.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.second.rankDisp() + << " " << tracks.second.zPosition() << " " << tracks.second.kSlope() << " " + << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; } } else if (patterns) { - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; } if (i < Nstubs2) { l1t::MuonStubRefVector stubs_proc; @@ -227,22 +237,25 @@ std::pair, std::vector > KMTF::proce if (tracks.first.id() & 0x1) edm::LogInfo("KMTF") << "1 " << (tracks.first.curvatureAtVertex() < 0 ? 1 : 0) << " " << tracks.first.ptPrompt() << " " << tracks.first.phiAtMuon() / (1 << 5) << " " - << tracks.first.coarseEta() << " " << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.first.rankPrompt() << " " << std::flush; + << int(round(tracks.first.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.first.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.first.rankPrompt() + << " " << tracks.first.zPosition() << " " << tracks.first.kSlope() << " " << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; if (tracks.second.id() & 0x2) - edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtVertex() < 0 ? 1 : 0) << " " + edm::LogInfo("KMTF") << "1 " << (tracks.second.curvatureAtMuon() < 0 ? 1 : 0) << " " << tracks.second.ptDisplaced() << " " << tracks.second.phiAtMuon() / (1 << 5) << " " - << tracks.second.coarseEta() << " " << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) - << " " << tracks.second.rankDisp(); + << int(round(tracks.second.eta() / (M_PI / (1 << 12)))) << " " + << int(tracks.second.dxy() * ap_ufixed<8, 1>(1.606)) << " " << tracks.second.rankDisp() + << " " << tracks.second.zPosition() << " " << tracks.second.kSlope() << " " + << std::flush; else - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0" << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; } } else if (patterns) { - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 " << std::flush; - edm::LogInfo("KMTF") << "0 0 0 0 0 0 0"; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0 " << std::flush; + edm::LogInfo("KMTF") << "0 0 0 0 0 0 0 0 0"; } //Now the shift register emulation in C_++ if (stubs4.size() > 1) { @@ -301,6 +314,8 @@ void KMTF::overlapCleanTrack(l1t::KMTFTrack& source, const l1t::KMTFTrack& other int id2 = vertex ? other.id() & 0x1 : other.id() & 0x2; bool keep = true; unsigned int pattern = 0; + unsigned int thetaPattern = 0; + if (id1 == 0) keep = false; else if (id1 != 0 && id2 != 0) { @@ -320,10 +335,13 @@ void KMTF::overlapCleanTrack(l1t::KMTFTrack& source, const l1t::KMTFTrack& other if (ok) { stubs.push_back(s1); pattern = pattern | (1 << (s1->depthRegion() - 1)); + if (s1->depthRegion() != 4 && s1->etaQuality() > 0) + thetaPattern |= (1 << (s1->depthRegion() - 1)); } } source.setStubs(stubs); source.setHitPattern(pattern); + source.setThetaDigiPattern(thetaPattern); } std::vector KMTF::cleanRegion(const std::vector& tracks2, diff --git a/L1Trigger/Phase2L1GMT/src/KMTFCore.cc b/L1Trigger/Phase2L1GMT/src/KMTFCore.cc index fcbca8ba7284c..0cb775261d990 100644 --- a/L1Trigger/Phase2L1GMT/src/KMTFCore.cc +++ b/L1Trigger/Phase2L1GMT/src/KMTFCore.cc @@ -1,7 +1,8 @@ #include "L1Trigger/Phase2L1GMT/interface/KMTFCore.h" using namespace Phase2L1GMT; KMTFCore::KMTFCore(const edm::ParameterSet& settings) - : lutService_(new KMTFLUTs(settings.getParameter("lutFile"))), + : lutService_(new KMTFLUTs(settings.getParameter("lutFile"), + settings.getParameter("lutThetaFile"))), verbose_(settings.getParameter("verbose")), initK_(settings.getParameter >("initialK")), initK2_(settings.getParameter >("initialK2")), @@ -47,13 +48,18 @@ KMTFCore::KMTFCore(const edm::ParameterSet& settings) combos1_(settings.getParameter >("combos1")), useOfflineAlgo_(settings.getParameter("useOfflineAlgo")), + zdeltaR_(settings.getParameter >("zdeltaR")), mScatteringPhi_(settings.getParameter >("mScatteringPhi")), mScatteringPhiB_(settings.getParameter >("mScatteringPhiB")), + mScatteringz_(settings.getParameter >("mScatteringz")), + mScatteringkSlope_(settings.getParameter >("mScatteringkSlope")), pointResolutionPhi_(settings.getParameter("pointResolutionPhi")), pointResolutionPhiB_(settings.getParameter("pointResolutionPhiB")), pointResolutionPhiBH_(settings.getParameter >("pointResolutionPhiBH")), pointResolutionPhiBL_(settings.getParameter >("pointResolutionPhiBL")), pointResolutionVertex_(settings.getParameter("pointResolutionVertex")), + pointResolutionz_(settings.getParameter >("pointResolutionz")), + pointResolutionkSlope_(settings.getParameter >("pointResolutionkSlope")), curvResolution1_(settings.getParameter >("curvResolution1")), curvResolution2_(settings.getParameter >("curvResolution2")) {} @@ -78,11 +84,10 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef default: throw cms::Exception("KMTFCore") << "Something really bad happend\n"; } - - l1t::KMTFTrack nullTrack(seed, seed->coord1(), correctedPhiB(seed)); + l1t::KMTFTrack nullTrack(seed, seed->coord1(), correctedPhiB(seed), seed->eta1(), satKSlope(seed->eta2())); seedQual = seed->quality(); for (const auto& mask : combinatorics) { - l1t::KMTFTrack track(seed, seed->coord1(), correctedPhiB(seed)); + l1t::KMTFTrack track(seed, seed->coord1(), correctedPhiB(seed), seed->eta1(), satKSlope(seed->eta2())); int phiB = correctedPhiB(seed); int charge; if (phiB == 0) @@ -103,16 +108,20 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef initialK = pow(2, BITSCURV - 1) - 1; if (initialK <= -pow(2, BITSCURV - 1)) initialK = -pow(2, BITSCURV - 1) + 1; - track.setCoordinates(seed->depthRegion(), initialK, seed->coord1(), phiB); + + track.setCoordinates(seed->depthRegion(), initialK, seed->coord1(), phiB, seed->eta1(), satKSlope(seed->eta2())); + if (seed->quality() < 6) { - track.setCoordinates(seed->depthRegion(), initialK, seed->coord1(), 0); + track.setCoordinates(seed->depthRegion(), initialK, seed->coord1(), 0, seed->eta1(), satKSlope(seed->eta2())); } if (verbose_) { - edm::LogInfo("KMTFCore") << "Initial state: phiB=" << phiB << " addr=" << address << " K=" << initialK; + edm::LogInfo("KMTFCore") << "Initial state: phiB=" << phiB << " addr=" << address << " K=" << initialK + << " z=" << seed->eta1() << " kSlope=" << satKSlope(seed->eta2()); } track.setHitPattern(hitPattern(track)); + track.setThetaDigiPattern(thetaDigiPattern(track)); //set covariance - l1t::KMTFTrack::CovarianceMatrix covariance; + l1t::CovarianceMatrix5dim covariance; float DK = curvResolution1_[track.step() - 1] + curvResolution2_[track.step() - 1] * initialK * initialK; if (seed->quality() < 6) DK = pow(2, 22); @@ -120,11 +129,25 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef covariance(0, 0) = DK * 4; covariance(0, 1) = 0; covariance(0, 2) = 0; + covariance(0, 3) = 0; + covariance(0, 4) = 0; covariance(1, 0) = 0; covariance(1, 1) = float(pointResolutionPhi_); covariance(1, 2) = 0; + covariance(1, 3) = 0; + covariance(1, 4) = 0; covariance(2, 0) = 0; covariance(2, 1) = 0; + covariance(2, 3) = 0; + covariance(2, 4) = 0; + covariance(3, 0) = 0; + covariance(3, 1) = 0; + covariance(3, 2) = 0; + covariance(3, 4) = 0; + covariance(4, 0) = 0; + covariance(4, 1) = 0; + covariance(4, 2) = 0; + covariance(4, 3) = 0; if (!(mask == 1 || mask == 2 || mask == 3 || mask == 4 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)) covariance(2, 2) = float(pointResolutionPhiB_); @@ -134,12 +157,19 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef else covariance(2, 2) = float(pointResolutionPhiBH_[seed->depthRegion() - 1]); } + if (seed->etaQuality() <= 0) { + covariance(3, 3) = pow(2, 22); + covariance(4, 4) = pow(2, 22); + } else { + covariance(3, 3) = float(pointResolutionz_[seed->depthRegion() - 1]); + covariance(4, 4) = float(pointResolutionkSlope_[seed->depthRegion() - 1]); + } track.setCovariance(covariance); - // if (verbose_) { edm::LogInfo("KMTFCore") << "New Kalman fit staring at step=" << track.step() << ", phi=" << track.positionAngle() - << ", phiB=" << track.bendingAngle() << ", with curvature=" << track.curvature(); + << ", phiB=" << track.bendingAngle() << ", z=" << track.zPosition() + << ",kSlope = " << track.kSlope() << ", with curvature=" << track.curvature(); edm::LogInfo("KMTFCore") << "BITMASK:" << std::flush; for (unsigned int i = 0; i < 4; ++i) edm::LogInfo("KMTFCore") << getBit(mask, i) << std::flush; @@ -152,8 +182,8 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef edm::LogInfo("KMTFCore") << "station=" << stub->depthRegion() << " phi=" << stub->coord1() << " phiB=" << correctedPhiB(stub) << " qual=" << stub->quality() << " tag=" << stub->id() << " sector=" << stub->phiRegion() - << " wheel=" << stub->etaRegion() << " fineEta= " << stub->eta1() << " " - << stub->eta2(); + << " wheel=" << stub->etaRegion() << " z= " << stub->eta1() + << " kSlope = " << satKSlope(stub->eta2()); edm::LogInfo("KMTFCore") << "------------------------------------------------------"; edm::LogInfo("KMTFCore") << "------------------------------------------------------"; } @@ -163,7 +193,8 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef while (track.step() > 0) { // muon station 1 if (track.step() == 1) { - track.setCoordinatesAtMuon(track.curvature(), track.positionAngle(), track.bendingAngle()); + track.setCoordinatesAtMuon( + track.curvature(), track.positionAngle(), track.bendingAngle(), track.zPosition(), track.kSlope()); passedU = estimateChiSquare(track, false); setRank(track, false); @@ -174,36 +205,64 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef setFourVectors(track); //calculate coarse eta ////////////////////// - if (verbose_) + if (verbose_) { edm::LogInfo("KMTFCore") << "Unconstrained PT in Muon System: pt=" << track.displacedP4().pt(); + } calculateEta(track); } propagate(track); if (verbose_) edm::LogInfo("KMTFCore") << "propagated Coordinates step:" << track.step() << "phi=" << track.positionAngle() - << "phiB=" << track.bendingAngle() << "K=" << track.curvature(); - if (track.step() > 0) + << "phiB=" << track.bendingAngle() << "K=" << track.curvature() + << ", z=" << track.zPosition() << ", kSlope=" << track.kSlope(); + if (track.step() > 0) { if (getBit(mask, track.step() - 1)) { std::pair bestStub = match(seed, stubs, track.step()); if (verbose_) edm::LogInfo("KMTFCore") << "Found match =" << bestStub.first << " index=" << bestStub.second << " number of all stubs=" << stubs.size(); - - if ((!bestStub.first) || (!update(track, stubs[bestStub.second], mask, seedQual))) + const bool updated = bestStub.first ? update(track, stubs[bestStub.second], mask, seedQual) : false; + if ((!bestStub.first) || (!updated)) { + if (verbose_) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG] Breaking Kalman chain before final ID assignment" + << " seed(st=" << seed->depthRegion() << ",sec=" << seed->phiRegion() + << ",wh=" << seed->etaRegion() << ",bx=" << seed->bxNum() << ",id=" << seed->id() + << ")" + << " step=" << track.step() << " mask=" << mask + << " hitPattern=" << track.hitPattern() << " matchFound=" << bestStub.first + << " matchIndex=" << bestStub.second << " updateOK=" << updated + << " partialID=" << track.id() << " nStubs=" << track.stubs().size(); + if (seed->depthRegion() == 2) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG_STA2] break-on-no-update/no-match" + << " seed(st=2,sec=" << seed->phiRegion() << ",wh=" << seed->etaRegion() + << ",bx=" << seed->bxNum() << ",id=" << seed->id() << ")" + << " currentStep=" << track.step() << " mask=" << mask + << " hitPattern=" << track.hitPattern() + << " thetaPattern=" << track.thetaDigiPattern() << " passU_sofar=" << passedU + << " passV_sofar=" << passedV << " id_sofar=" << track.id() + << " phiMuon=" << track.phiAtMuon() << " phiVertex=" << track.phiAtVertex() + << " z=" << track.zPosition() << " kSlope=" << track.kSlope(); + } + } break; + } if (verbose_) { edm::LogInfo("KMTFCore") << "updated Coordinates step:" << track.step() << " phi=" << track.positionAngle() - << " phiB=" << track.bendingAngle() << " K=" << track.curvature(); + << " phiB=" << track.bendingAngle() << " K=" << track.curvature() + << ", z=" << track.zPosition() << ", kSlope=" << track.kSlope(); } } + } if (track.step() == 0) { - track.setCoordinatesAtVertex(track.curvature(), track.positionAngle(), track.bendingAngle()); + track.setCoordinatesAtVertex( + track.curvature(), track.positionAngle(), track.bendingAngle(), track.zPosition(), track.kSlope()); if (verbose_) edm::LogInfo("KMTFCore") << " Coordinates before vertex constraint step:" << track.step() << " phi=" << track.phiAtVertex() << " dxy=" << track.dxy() - << " K=" << track.curvatureAtVertex(); + << " K=" << track.curvatureAtVertex() << ", z=" << track.zPosition() + << ", kSlope=" << track.kSlope(); //apply vertex constraint for non single tracks if (track.stubs().size() > 1) @@ -226,12 +285,37 @@ std::pair KMTFCore::chain(const l1t::MuonStubRef } setFourVectors(track); //finally set the displaced or prompt ID + const int oldId = track.id(); track.setIDFlag(passedV, passedU); + if (verbose_) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG] setIDFlag transition" + << " seed(st=" << seed->depthRegion() << ",sec=" << seed->phiRegion() + << ",wh=" << seed->etaRegion() << ",bx=" << seed->bxNum() << ",id=" << seed->id() + << ")" + << " mask=" << mask << " hitPattern=" << track.hitPattern() + << " thetaPattern=" << track.thetaDigiPattern() << " passPrompt=" << passedV + << " passDisp=" << passedU << " promptChi2=" << track.approxPromptChi2() + << " dispChi2=" << track.approxDispChi2() << " rankPrompt=" << track.rankPrompt() + << " rankDisp=" << track.rankDisp() << " oldID=" << oldId << " newID=" << track.id() + << " phiMuon=" << track.phiAtMuon() << " phiVertex=" << track.phiAtVertex() + << " nStubs=" << track.stubs().size(); + if (seed->depthRegion() == 2) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG_STA2] final-id-decision" + << " seed(st=2,sec=" << seed->phiRegion() << ",wh=" << seed->etaRegion() + << ",bx=" << seed->bxNum() << ",id=" << seed->id() << ")" + << " mask=" << mask << " passPrompt=" << passedV << " passDisp=" << passedU + << " oldID=" << oldId << " newID=" << track.id() + << " rankPrompt=" << track.rankPrompt() << " rankDisp=" << track.rankDisp() + << " promptChi2=" << track.approxPromptChi2() + << " dispChi2=" << track.approxDispChi2() << " phiMuon=" << track.phiAtMuon() + << " phiVertex=" << track.phiAtVertex() << " z=" << track.zPosition() + << " kSlope=" << track.kSlope() << " nStubs=" << track.stubs().size(); + } + } if (verbose_) edm::LogInfo("KMTFCore") << "Floating point coordinates at vertex: pt=" << track.pt() << ", eta=" << track.eta() << " phi=" << track.phi(); - pretracks.push_back(track); } } @@ -281,8 +365,15 @@ std::vector KMTFCore::clean(const std::vector& t for (const auto& track : tracks) { if (vertex) { - if ((track.id() & 0x1) == 0) + if ((track.id() & 0x1) == 0) { + if (verbose_) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG] clean(prompt) dropping pretrack due to prompt ID bit=0" + << " seedStep=" << seed << " id=" << track.id() + << " hitPattern=" << track.hitPattern() << " rankPrompt=" << track.rankPrompt() + << " rankDisp=" << track.rankDisp() << " nStubs=" << track.stubs().size(); + } continue; + } if (verbose_) edm::LogInfo("KMTFCore") << "Chain Cleaning : Pre Track = pattern = " << track.rankPrompt() << " rank=" << track.hitPattern(); @@ -290,8 +381,15 @@ std::vector KMTFCore::clean(const std::vector& t infoTrack[track.hitPattern()] = track; } else { - if ((track.id() & 0x2) == 0) + if ((track.id() & 0x2) == 0) { + if (verbose_) { + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG] clean(displaced) dropping pretrack due to displaced ID bit=0" + << " seedStep=" << seed << " id=" << track.id() + << " hitPattern=" << track.hitPattern() << " rankPrompt=" << track.rankPrompt() + << " rankDisp=" << track.rankDisp() << " nStubs=" << track.stubs().size(); + } continue; + } infoRank[track.hitPattern()] = track.rankDisp(); infoTrack[track.hitPattern()] = track; } @@ -331,9 +429,21 @@ std::vector KMTFCore::clean(const std::vector& t selected = 1; auto search = infoTrack.find(selected); - if (search != infoTrack.end()) + if (search != infoTrack.end()) { + if (verbose_) { + const auto& kept = search->second; + edm::LogInfo("KMTFCore") << "[KMTF_VALIDDBG] clean(" << (vertex ? "prompt" : "displaced") + << ") keeping selected pretrack" + << " seedStep=" << seed << " selectedPattern=" << selected << " id=" << kept.id() + << " rankPrompt=" << kept.rankPrompt() << " rankDisp=" << kept.rankDisp() + << " hitPattern=" << kept.hitPattern() << " thetaPattern=" << kept.thetaDigiPattern() + << " phiMuon=" << kept.phiAtMuon() << " phiVertex=" << kept.phiAtVertex() + << " curvMuon=" << kept.curvatureAtMuon() << " curvVertex=" << kept.curvatureAtVertex() + << " z=" << kept.zPosition() << " kSlope=" << kept.kSlope() + << " nStubs=" << kept.stubs().size(); + } out.push_back(search->second); - + } return out; } @@ -429,8 +539,10 @@ void KMTFCore::propagate(l1t::KMTFTrack& track) { int phi = track.positionAngle(); int phiB = track.bendingAngle(); unsigned int step = track.step(); - + int z = track.zPosition(); + int kSlope = track.kSlope(); int charge = 1; + if (K != 0) charge = K / fabs(K); @@ -477,6 +589,27 @@ void KMTFCore::propagate(l1t::KMTFTrack& track) { << ap_ufixed(bPhiB_[step - 1]).to_float() << " = " << phiB12.to_int(); } + //z and kSlope propagation + // dR conversion factor comes from ZRES_CONV and KRES_CONV: ZRES_CONV/KRES_CONV=(65536/1500)/(65536/2) to go from physical (cm) to digitized units + double zdR_CONV = 2.0 / 1500.0; + double zdeltaR_dig = zdeltaR_[step - 1] * zdR_CONV; + int kSlopeNew = kSlope; + int zNew; + zNew = + ap_fixed(ap_fixed(z) - ap_ufixed(zdeltaR_dig) * + ap_fixed(kSlope)); + if (zNew > (1 << (BITSZ - 1)) - 1) { + if (verbose_) { + edm::LogInfo("KMTFCore") << "z saturated high during propagation, step=" << step; + } + zNew = (1 << (BITSZ - 1)) - 1; + } else if (zNew < -(1 << (BITSZ - 1))) { + if (verbose_) { + edm::LogInfo("KMTFCore") << "z saturated low during propagation, step=" << step; + } + zNew = -(1 << (BITSZ - 1)); + } + //Only for the propagation to vertex we use second order; if (step == 1) { ap_fixed<10, 4> aPhiB = aPhiB_[step - 1]; @@ -504,51 +637,76 @@ void KMTFCore::propagate(l1t::KMTFTrack& track) { //where we want to check what is happening in the covariance matrix //Create the transformation matrix - double a[9]; - a[0] = 1.; + double a[25]; + a[0] = 1.0; a[1] = 0.0; a[2] = 0.0; - a[3] = aPhi_[step - 1]; - // a[3] = 0.0; - a[4] = 1.0; - a[5] = -bPhi_[step - 1]; - //a[6]=0.0; - a[6] = aPhiB_[step - 1]; + a[3] = 0.0; + a[4] = 0.0; + a[5] = aPhi_[step - 1]; + a[6] = 1.0; + a[7] = -bPhi_[step - 1]; + a[8] = 0.0; + a[9] = 0.0; + a[10] = aPhiB_[step - 1]; if (step == 1) - a[6] = aPhiB_[step - 1] / 2.0; - - a[7] = 0.0; - a[8] = bPhiB_[step - 1]; - - ROOT::Math::SMatrix P(a, 9); + a[10] = aPhiB_[step - 1] / 2.0; + a[11] = 0.0; + a[12] = bPhiB_[step - 1]; + a[13] = 0.0; + a[14] = 0.0; + a[15] = 0.0; + a[16] = 0.0; + a[17] = 0.0; + a[18] = 1.0; + a[19] = -zdeltaR_dig; + a[20] = 0.0; + a[21] = 0.0; + a[22] = 0.0; + a[23] = 0.0; + a[24] = 1.0; + + ROOT::Math::SMatrix P(a, 25); const std::vector& covLine = track.covariance(); - l1t::KMTFTrack::CovarianceMatrix cov(covLine.begin(), covLine.end()); + l1t::CovarianceMatrix5dim cov(covLine.begin(), covLine.end()); cov = ROOT::Math::Similarity(P, cov); //Add the multiple scattering double phiRMS = mScatteringPhi_[step - 1] * K * K; double phiBRMS = mScatteringPhiB_[step - 1] * K * K; + double zRMS = mScatteringz_[step - 1] * K * K; + double kSlopeRMS = mScatteringkSlope_[step - 1] * K * K; - std::vector b(6); + std::vector b(15); b[0] = 0; b[1] = 0; b[2] = phiRMS; b[3] = 0; b[4] = 0; b[5] = phiBRMS; - - reco::Candidate::CovarianceMatrix MS(b.begin(), b.end()); + b[6] = 0; + b[7] = 0; + b[8] = 0; + b[9] = zRMS; + b[10] = 0; + b[11] = 0; + b[12] = 0; + b[13] = 0; + b[14] = kSlopeRMS; + + l1t::CovarianceMatrix5dim MS(b.begin(), b.end()); cov = cov + MS; if (verbose_) { edm::LogInfo("KMTFCore") << "Covariance term for phiB = " << cov(2, 2); edm::LogInfo("KMTFCore") << "Multiple scattering term for phiB = " << MS(2, 2); + edm::LogInfo("KMTFCore") << "Multiple scattering term for z = " << MS(3, 3); + edm::LogInfo("KMTFCore") << "Multiple scattering term for kSlope = " << MS(4, 4); } - track.setCovariance(cov); - track.setCoordinates(step - 1, KNew, phiNew, phiBNew); + track.setCoordinates(step - 1, KNew, phiNew, phiBNew, zNew, kSlopeNew); } bool KMTFCore::update(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, int mask, int seedQual) { @@ -567,77 +725,151 @@ bool KMTFCore::updateOffline(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub int trackK = track.curvature(); int trackPhi = track.positionAngle(); int trackPhiB = track.bendingAngle(); + int trackz = track.zPosition(); + int trackSlope = track.kSlope(); int phi = stub->coord1(); int phiB = correctedPhiB(stub); + int z = stub->eta1(); + int kSlope = satKSlope(stub->eta2()); + int priorThetaPattern = track.thetaDigiPattern(); + int priorPhiPattern = track.hitPattern(); + int seedStation = track.stubs().empty() ? 0 : track.stubs()[0]->depthRegion(); - Vector2 residual; + Vector4 residual; residual[0] = ap_fixed(phi - trackPhi); residual[1] = phiB - trackPhiB; + residual[2] = ap_fixed(z - trackz); + residual[3] = ap_fixed(kSlope - trackSlope); - Matrix23 H; + Matrix45 H; H(0, 0) = 0.0; H(0, 1) = 1.0; H(0, 2) = 0.0; + H(0, 3) = 0.0; + H(0, 4) = 0.0; H(1, 0) = 0.0; H(1, 1) = 0.0; H(1, 2) = 1.0; + H(1, 3) = 0.0; + H(1, 4) = 0.0; + H(2, 0) = 0.0; + H(2, 1) = 0.0; + H(2, 2) = 0.0; + H(2, 3) = 1.0; + H(2, 4) = 0.0; + H(3, 0) = 0.0; + H(3, 1) = 0.0; + H(3, 2) = 0.0; + H(3, 3) = 0.0; + H(3, 4) = 1.0; const auto r11{stub->quality() < 6 ? pointResolutionPhiBL_[track.step() - 1] : pointResolutionPhiBH_[track.step() - 1]}; - const double r[]{pointResolutionPhi_, 0.0, 0.0, r11}; - const CovarianceMatrix2 R(r, 4); - + const double r[]{pointResolutionPhi_, + 0.0, + r11, + 0.0, + 0.0, + pointResolutionz_[track.step() - 1], + 0.0, + 0.0, + 0.0, + pointResolutionkSlope_[track.step() - 1]}; + + const CovarianceMatrix4 R(r, 10); const std::vector& covLine = track.covariance(); - const l1t::KMTFTrack::CovarianceMatrix cov(covLine.begin(), covLine.end()); - CovarianceMatrix2 S = ROOT::Math::Similarity(H, cov) + R; + const l1t::CovarianceMatrix5dim cov(covLine.begin(), covLine.end()); + + CovarianceMatrix4 S = ROOT::Math::Similarity(H, cov) + R; if (!S.Invert()) return false; - const Matrix32 Gain = cov * ROOT::Math::Transpose(H) * S; + const Matrix54 Gain = cov * ROOT::Math::Transpose(H) * S; - track.setKalmanGain( - track.step(), fabs(trackK), Gain(0, 0), Gain(0, 1), Gain(1, 0), Gain(1, 1), Gain(2, 0), Gain(2, 1)); + //track.setKalmanGain( + //track.step(), fabs(trackK), Gain(0, 0), Gain(0, 1), 1, 0, Gain(2, 0), Gain(2, 1)); int KNew = (trackK + int(Gain(0, 0) * residual(0) + Gain(0, 1) * residual(1))); if (fabs(KNew) > pow(2, BITSCURV - 1)) return false; - int phiNew = wrapAround(trackPhi + residual(0), pow(2, BITSPHI - 1)); int phiBNew = wrapAround(trackPhiB + int(Gain(2, 0) * residual(0) + Gain(2, 1) * residual(1)), pow(2, BITSPHIB - 1)); + int zNew = trackz + int(Gain(3, 2) * residual(2) + Gain(3, 3) * residual(3)); + int kSlopeNew = trackSlope + int(Gain(4, 2) * residual(2) + Gain(4, 3) * residual(3)); - track.setResidual(stub->depthRegion() - 1, fabs(phi - phiNew) + fabs(phiB - phiBNew)); + if ((zNew > ((1 << (BITSZ - 1)) - 1)) || (zNew < (-(1 << (BITSZ - 1))))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "z saturated in updateOffline"; + return false; + } + if ((kSlopeNew > ((1 << (BITSKSLOPE - 1)) - 1)) || (kSlopeNew < (-(1 << (BITSKSLOPE - 1))))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "kSlope saturated in updateOffline"; + return false; + } + + track.setResidual(stub->depthRegion() - 1, + fabs(phi - phiNew) + fabs(phiB - phiBNew) + fabs(z - zNew) + fabs(kSlope - kSlopeNew)); if (verbose_) { - edm::LogInfo("KMTFCore") << "residual " << phi << " - " << trackPhi << " = " << int(residual[0]) << " " << phiB - << " - " << trackPhiB << " = " << int(residual[1]); - edm::LogInfo("KMTFCore") << "Gains offline: " << Gain(0, 0) << " " << Gain(0, 1) << " " << Gain(2, 0) << " " - << Gain(2, 1); - edm::LogInfo("KMTFCore") << " K = " << trackK << " + " << Gain(0, 0) << " * " << residual(0) << " + " << Gain(0, 1) - << " * " << residual(1); - edm::LogInfo("KMTFCore") << " phiB = " << trackPhiB << " + " << Gain(2, 0) << " * " << residual(0) << " + " - << Gain(2, 1) << " * " << residual(1); + edm::LogInfo("KMTFCore") << "residual(0): " << phi << "-" << trackPhi << " = " << residual(0) << "\n"; + edm::LogInfo("KMTFCore") << "residual(1): " << phiB << "-" << trackPhiB << " = " << residual(1) << "\n"; + edm::LogInfo("KMTFCore") << "residual(2): " << z << "-" << trackz << " = " << residual(2) << "\n"; + edm::LogInfo("KMTFCore") << "residual(3): " << kSlope << "-" << trackSlope << " = " << residual(3) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(0,0): " << Gain(0, 0) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(0,1): " << Gain(0, 1) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(2,0): " << Gain(2, 0) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(2,1): " << Gain(2, 1) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(3,2): " << Gain(3, 2) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(3,3): " << Gain(3, 3) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(4,2): " << Gain(4, 2) << "\n"; + edm::LogInfo("KMTFCore") << "Gain(4,3): " << Gain(4, 3) << "\n"; + edm::LogInfo("KMTFCore") << " KNew = " << trackK << "+" << Gain(0, 0) << "*" << residual(0) << "+" << Gain(0, 1) + << "*" << residual(1) << " = " << KNew << "\n"; + edm::LogInfo("KMTFCore") << " phiNew = " << trackPhi << "+" << residual(0) << " = " << phiNew << "\n"; + edm::LogInfo("KMTFCore") << " phiBNew = " << trackPhiB << "+" << Gain(2, 0) << "*" << residual(0) << "+" + << Gain(2, 1) << "*" << residual(1) << " = " << phiBNew << "\n"; + edm::LogInfo("KMTFCore") << " zNew = " << trackz << "+" << Gain(3, 2) << "*" << residual(2) << "+" << Gain(3, 3) + << "*" << residual(3) << " = " << zNew << "\n"; + edm::LogInfo("KMTFCore") << " kSlopeNew = " << trackSlope << "+" << Gain(4, 2) << "*" << residual(2) << "+" + << Gain(4, 3) << "*" << residual(3) << " = " << kSlopeNew << "\n"; } - track.setCoordinates(track.step(), KNew, phiNew, phiBNew); - const auto covNew{cov - Gain * (H * cov)}; - l1t::KMTFTrack::CovarianceMatrix c; + track.setCoordinates(track.step(), KNew, phiNew, phiBNew, zNew, kSlopeNew); + const Matrix55 covNew = cov - Gain * (H * cov); + l1t::CovarianceMatrix5dim c; - c(0, 0) = covNew(0, 0); - c(0, 1) = covNew(0, 1); - c(0, 2) = covNew(0, 2); - c(1, 0) = covNew(1, 0); - c(1, 1) = covNew(1, 1); - c(1, 2) = covNew(1, 2); - c(2, 0) = covNew(2, 0); - c(2, 1) = covNew(2, 1); - c(2, 2) = covNew(2, 2); - if (verbose_) { - edm::LogInfo("KMTFCore") << "Post Fit Covariance Matrix " << cov(0, 0) << " " << cov(1, 1) << " " << cov(2, 2); + for (int i = 0; i < 5; i++) { + for (int j = i; j < 5; j++) { + c(i, j) = covNew(i, j); + } } track.setCovariance(c); track.addStub(stub); track.setHitPattern(hitPattern(track)); + track.setThetaDigiPattern(thetaDigiPattern(track)); + track.setThetaGain(track.step(), + fabs(trackK), + priorThetaPattern, + seedStation, + priorPhiPattern, + Gain(3, 2), + Gain(3, 3), + Gain(4, 2), + Gain(4, 3)); + + if (verbose_) { + edm::LogInfo("KMTFCore") << "step: " << track.step() << "\n"; + edm::LogInfo("KMTFCore") << "|K|: " << fabs(trackK) << "\n"; + edm::LogInfo("KMTFCore") << "priorThetaPattern: " << priorThetaPattern << "\n"; + edm::LogInfo("KMTFCore") << "priorPhiPattern: " << priorPhiPattern << "\n"; + edm::LogInfo("KMTFCore") << "seed Station: " << seedStation << "\n"; + edm::LogInfo("KMTFCore") << "Gain32: " << Gain(3, 2) << "\n"; + edm::LogInfo("KMTFCore") << "Gain33: " << Gain(3, 3) << "\n"; + edm::LogInfo("KMTFCore") << "Gain42: " << Gain(4, 2) << "\n"; + edm::LogInfo("KMTFCore") << "Gain43: " << Gain(4, 3) << "\n"; + } return true; } @@ -646,54 +878,106 @@ bool KMTFCore::updateOffline1D(l1t::KMTFTrack& track, const l1t::MuonStubRef& st int trackK = track.curvature(); int trackPhi = track.positionAngle(); int trackPhiB = track.bendingAngle(); - + int trackz = track.zPosition(); + int trackSlope = track.kSlope(); int phi = stub->coord1(); + int z = stub->eta1(); + int kSlope = satKSlope(stub->eta2()); + int priorThetaPattern = track.thetaDigiPattern(); + int priorPhiPattern = track.hitPattern(); + int seedStation = track.stubs().empty() ? 0 : track.stubs()[0]->depthRegion(); - double residual = ap_fixed(phi - trackPhi); + Vector3 residual; + residual[0] = ap_fixed(phi - trackPhi); + residual[1] = ap_fixed(z - trackz); + residual[2] = ap_fixed(kSlope - trackSlope); - if (verbose_) - edm::LogInfo("KMTFCore") << "residuals " << phi << " - " << trackPhi << " = " << int(residual); + if (verbose_) { + edm::LogInfo("KMTFCore") << "residual phi: " << phi << " - " << trackPhi << " = " << int(residual(0)); + edm::LogInfo("KMTFCore") << "residual z: " << z << " - " << trackz << " = " << int(residual(1)); + edm::LogInfo("KMTFCore") << "residual kSlope: " << kSlope << " - " << trackSlope << " = " << int(residual(2)); + } - Matrix13 H; + Matrix35 H; H(0, 0) = 0.0; H(0, 1) = 1.0; H(0, 2) = 0.0; + H(0, 3) = 0.0; + H(0, 4) = 0.0; + H(1, 0) = 0.0; + H(1, 1) = 0.0; + H(1, 2) = 0.0; + H(1, 3) = 1.0; + H(1, 4) = 0.0; + H(2, 0) = 0.0; + H(2, 1) = 0.0; + H(2, 2) = 0.0; + H(2, 3) = 0.0; + H(2, 4) = 1.0; + + const double r[]{ + pointResolutionPhi_, 0.0, pointResolutionz_[track.step() - 1], 0.0, 0.0, pointResolutionkSlope_[track.step() - 1]}; + const CovarianceMatrix3 R(r, 6); const std::vector& covLine = track.covariance(); - l1t::KMTFTrack::CovarianceMatrix cov(covLine.begin(), covLine.end()); + l1t::CovarianceMatrix5dim cov(covLine.begin(), covLine.end()); - double S = ROOT::Math::Similarity(H, cov)(0, 0) + pointResolutionPhi_; - if (S == 0.0) + CovarianceMatrix3 S = ROOT::Math::Similarity(H, cov) + R; + if (!S.Invert()) return false; - Matrix31 Gain = cov * ROOT::Math::Transpose(H) / S; + Matrix53 Gain = cov * ROOT::Math::Transpose(H) * S; - track.setKalmanGain(track.step(), fabs(trackK), Gain(0, 0), 0.0, Gain(1, 0), 0.0, Gain(2, 0), 0.0); + //track.setKalmanGain(track.step(), fabs(trackK), Gain(0, 0), 0.0, 1, 0, Gain(2, 0), 0.0); - int KNew = wrapAround(trackK + int(Gain(0, 0) * residual), pow(2, BITSCURV - 1)); - int phiNew = wrapAround(trackPhi + residual, pow(2, BITSPHI - 1)); - int phiBNew = wrapAround(trackPhiB + int(Gain(2, 0) * residual), pow(2, BITSPHIB - 1)); - track.setCoordinates(track.step(), KNew, phiNew, phiBNew); - Matrix33 covNew = cov - Gain * (H * cov); - l1t::KMTFTrack::CovarianceMatrix c; + int KNew = wrapAround(trackK + int(Gain(0, 0) * residual(0)), pow(2, BITSCURV - 1)); + int phiNew = wrapAround(trackPhi + residual(0), pow(2, BITSPHI - 1)); + int phiBNew = wrapAround(trackPhiB + int(Gain(2, 0) * residual(0)), pow(2, BITSPHIB - 1)); + int zNew = trackz + int(Gain(3, 1) * residual(1) + Gain(3, 2) * residual(2)); + int kSlopeNew = trackSlope + int(Gain(4, 1) * residual(1) + Gain(4, 2) * residual(2)); + + if ((zNew > (1 << (BITSZ - 1)) - 1) || (zNew < -(pow(2, BITSZ - 1)))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "z saturated in updateOffline1D"; + return false; + } + if ((kSlopeNew > (pow(2, BITSKSLOPE - 1) - 1)) || (kSlopeNew < -(pow(2, BITSKSLOPE - 1)))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "kSlope saturated in updateOffline1D"; + return false; + } + + track.setCoordinates(track.step(), KNew, phiNew, phiBNew, zNew, kSlopeNew); + + Matrix55 covNew = cov - Gain * (H * cov); + l1t::CovarianceMatrix5dim c; if (verbose_) { - edm::LogInfo("KMTFCore") << "phiUpdate: " << int(Gain(0, 0) * residual) << " " << int(Gain(2, 0) * residual); - edm::LogInfo("KMTFCore") << " K = " << trackK << " + " << Gain(0, 0) << " * " << residual; - edm::LogInfo("KMTFCore") << " phiBNew = " << trackPhiB << " + " << Gain(2, 0) << " * " << residual; + edm::LogInfo("KMTFCore") << " phiNew: " << trackPhi << "+" << residual(0); + edm::LogInfo("KMTFCore") << " KNew = " << trackK << "+" << Gain(0, 0) << "*" << residual(0); + edm::LogInfo("KMTFCore") << " phiBNew = " << trackPhiB << "+" << Gain(2, 0) << "*" << residual(0); + edm::LogInfo("KMTFCore") << " zNew = " << trackz << "+" << Gain(3, 1) << "*" << residual(1) << "+" << Gain(3, 2) + << "*" << residual(2); + edm::LogInfo("KMTFCore") << " kSlopeNew = " << trackSlope << "+" << Gain(4, 1) << "*" << residual(1) << "+" + << Gain(4, 2) << "*" << residual(2); } - c(0, 0) = covNew(0, 0); - c(0, 1) = covNew(0, 1); - c(0, 2) = covNew(0, 2); - c(1, 0) = covNew(1, 0); - c(1, 1) = covNew(1, 1); - c(1, 2) = covNew(1, 2); - c(2, 0) = covNew(2, 0); - c(2, 1) = covNew(2, 1); - c(2, 2) = covNew(2, 2); + for (int i = 0; i < 5; i++) + for (int j = i; j < 5; j++) + c(i, j) = covNew(i, j); + track.setCovariance(c); track.addStub(stub); track.setHitPattern(hitPattern(track)); + track.setThetaDigiPattern(thetaDigiPattern(track)); + track.setThetaGain1D(track.step(), + fabs(trackK), + priorThetaPattern, + seedStation, + priorPhiPattern, + Gain(3, 1), + Gain(3, 2), + Gain(4, 1), + Gain(4, 2)); return true; } @@ -702,13 +986,19 @@ bool KMTFCore::updateLUT(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, in int trackK = track.curvature(); int trackPhi = track.positionAngle(); int trackPhiB = track.bendingAngle(); + int trackz = track.zPosition(); + int trackSlope = track.kSlope(); int phi = stub->coord1(); int phiB = correctedPhiB(stub); + int z = stub->eta1(); + int kSlope = satKSlope(stub->eta2()); - Vector2 residual; + Vector4 residual; ap_fixed residualPhi = phi - trackPhi; ap_fixed residualPhiB = phiB - trackPhiB; + ap_fixed residualz = z - trackz; + ap_fixed residualSlope = kSlope - trackSlope; if (verbose_) edm::LogInfo("KMTFCore") << "residual " << phi << " - " << trackPhi << " = " << residualPhi.to_int() << " " << phiB @@ -719,6 +1009,7 @@ bool KMTFCore::updateLUT(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, in absK = pow(2, BITSCURV - 2) - 1; std::vector GAIN; + std::vector GAIN_THETA; if (verbose_) { edm::LogInfo("KMTFCore") << "Looking up LUTs for mask=" << mask << " with hit pattern=" << track.hitPattern(); } @@ -731,6 +1022,14 @@ bool KMTFCore::updateLUT(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, in } else { GAIN = lutService_->trackGain2(track.step(), track.hitPattern(), absK / 32, seedQual, stub->quality()); } + + const bool is2D = (mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12); + if (!(is2D) && (track.step() == 2) && (track.thetaDigiPattern() == 0)) { + GAIN_THETA = lutService_->trackGainTheta2(track.step(), track.thetaDigiPattern(), track.hitPattern(), absK / 32); + } else { + GAIN_THETA = lutService_->trackGainTheta(track.step(), track.thetaDigiPattern(), absK / 32, is2D); + } + if (verbose_) { edm::LogInfo("KMTFCore") << "Gains (fp): " << GAIN[0] << " " << GAIN[1] << " " << GAIN[2] << " " << GAIN[3]; @@ -745,8 +1044,6 @@ bool KMTFCore::updateLUT(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, in << ap_ufixed(GAIN[3]).to_float(); } - track.setKalmanGain(track.step(), fabs(trackK), GAIN[0], GAIN[1], 1, 0, GAIN[2], GAIN[3]); - int KNew; if (!(mask == 3 || mask == 5 || mask == 9 || mask == 6 || mask == 10 || mask == 12)) { KNew = ap_fixed(ap_fixed(trackK) + @@ -794,9 +1091,31 @@ bool KMTFCore::updateLUT(l1t::KMTFTrack& track, const l1t::MuonStubRef& stub, in if ((phiBNew > (pow(2, BITSPHIB - 1) - 1)) || (phiBNew < (-pow(2, BITSPHIB - 1)))) return false; - track.setCoordinates(track.step(), KNew, phiNew, phiBNew); + int zNew; + int kSlopeNew; + zNew = ap_fixed(ap_fixed(trackz) + + ap_ufixed(GAIN_THETA[0]) * residualz - + ap_ufixed(GAIN_THETA[1]) * residualSlope); + kSlopeNew = ap_fixed(ap_fixed(trackSlope) - + ap_ufixed(GAIN_THETA[2]) * residualz + + ap_ufixed(GAIN_THETA[3]) * residualSlope); + + if ((zNew > (1 << (BITSZ - 1)) - 1) || (zNew < -(pow(2, BITSZ - 1)))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "z has saturated"; + return false; + } + if ((kSlopeNew > (pow(2, BITSKSLOPE - 1) - 1)) || (kSlopeNew < -(pow(2, BITSKSLOPE - 1)))) { + if (verbose_) + edm::LogInfo("KMTFCore") << "kSlope has saturated"; + return false; + } + + track.setCoordinates(track.step(), KNew, phiNew, phiBNew, zNew, kSlopeNew); track.addStub(stub); track.setHitPattern(hitPattern(track)); + track.setThetaDigiPattern(thetaDigiPattern(track)); + if (verbose_) { edm::LogInfo("KMTFCore") << "Stub station =" << stub->depthRegion(); @@ -815,17 +1134,19 @@ void KMTFCore::vertexConstraint(l1t::KMTFTrack& track) { void KMTFCore::vertexConstraintOffline(l1t::KMTFTrack& track) { double residual = -track.dxy(); - Matrix13 H; + Matrix15 H; H(0, 0) = 0; H(0, 1) = 0; H(0, 2) = 1; + H(0, 3) = 0; + H(0, 4) = 0; const std::vector& covLine = track.covariance(); - l1t::KMTFTrack::CovarianceMatrix cov(covLine.begin(), covLine.end()); + l1t::CovarianceMatrix5dim cov(covLine.begin(), covLine.end()); double S = (ROOT::Math::Similarity(H, cov))(0, 0) + pointResolutionVertex_; S = 1.0 / S; - Matrix31 Gain = cov * (ROOT::Math::Transpose(H)) * S; + Matrix51 Gain = cov * (ROOT::Math::Transpose(H)) * S; track.setKalmanGain(track.step(), fabs(track.curvature()), Gain(0, 0), Gain(1, 0), Gain(2, 0)); if (verbose_) { @@ -838,20 +1159,35 @@ void KMTFCore::vertexConstraintOffline(l1t::KMTFTrack& track) { int dxyNew = wrapAround(int(track.dxy() + Gain(2, 0) * residual), pow(2, BITSPHIB)); if (verbose_) edm::LogInfo("KMTFCore") << "Post fit impact parameter=" << dxyNew; - track.setCoordinatesAtVertex(KNew, phiNew, -residual); - Matrix33 covNew = cov - Gain * (H * cov); - l1t::KMTFTrack::CovarianceMatrix c; + track.setCoordinatesAtVertex(KNew, phiNew, -residual, track.zPosition(), track.kSlope()); + Matrix55 covNew = cov - Gain * (H * cov); + l1t::CovarianceMatrix5dim c; c(0, 0) = covNew(0, 0); c(0, 1) = covNew(0, 1); c(0, 2) = covNew(0, 2); + c(0, 3) = covNew(0, 3); + c(0, 4) = covNew(0, 4); c(1, 0) = covNew(1, 0); c(1, 1) = covNew(1, 1); c(1, 2) = covNew(1, 2); + c(1, 3) = covNew(1, 3); + c(1, 4) = covNew(1, 4); c(2, 0) = covNew(2, 0); c(2, 1) = covNew(2, 1); c(2, 2) = covNew(2, 2); + c(2, 3) = covNew(2, 3); + c(2, 4) = covNew(2, 4); + c(3, 0) = covNew(3, 0); + c(3, 1) = covNew(3, 1); + c(3, 2) = covNew(3, 2); + c(3, 3) = covNew(3, 3); + c(3, 4) = covNew(3, 4); + c(4, 0) = covNew(4, 0); + c(4, 1) = covNew(4, 1); + c(4, 2) = covNew(4, 2); + c(4, 3) = covNew(4, 3); + c(4, 4) = covNew(4, 4); track.setCovariance(c); - // track.covariance = track.covariance - Gain*H*track.covariance; } void KMTFCore::vertexConstraintLUT(l1t::KMTFTrack& track) { @@ -876,7 +1212,7 @@ void KMTFCore::vertexConstraintLUT(l1t::KMTFTrack& track) { //int p_0 = fp_product(GAIN.second, int(residual), 7); int p_0 = GAIN.second * int(residual); int phiNew = wrapAround(track.positionAngle() + p_0, pow(2, BITSPHI - 1)); - track.setCoordinatesAtVertex(KNew, phiNew, -residual); + track.setCoordinatesAtVertex(KNew, phiNew, -residual, track.zPosition(), track.kSlope()); } int KMTFCore::hitPattern(const l1t::KMTFTrack& track) { @@ -887,6 +1223,18 @@ int KMTFCore::hitPattern(const l1t::KMTFTrack& track) { return mask; } +//returns mask pattern based on stubs with a theta digi or not. e.g. 0010 is a track with stub at station 2 with a theta digi stub +int KMTFCore::thetaDigiPattern(const l1t::KMTFTrack& track) { + unsigned int mask = 0; + for (const auto& stub : track.stubs()) { + if (stub->depthRegion() == 4) + continue; + if (stub->etaQuality() > 0) + mask |= (1 << (stub->depthRegion() - 1)); + } + return mask; +} + int KMTFCore::customBitmask(unsigned int bit1, unsigned int bit2, unsigned int bit3, unsigned int bit4) { return bit1 * 1 + bit2 * 2 + bit3 * 4 + bit4 * 8; } @@ -894,7 +1242,15 @@ int KMTFCore::customBitmask(unsigned int bit1, unsigned int bit2, unsigned int b bool KMTFCore::getBit(int bitmask, int pos) { return (bitmask & (1 << pos)) >> pos; } void KMTFCore::setFourVectors(l1t::KMTFTrack& track) { - int etaINT = track.coarseEta(); + //int etaINT = track.coarseEta(); + //legacy version of KMTF used coarse Eta to set the eta in SAMuons based on wheel-based LUT + //new approach uses kSlope to build eta for SAMuons + //matched KMTF output to gen muons and performed linear fit of gen eta vs kSlope of KMTF output + //eta=m*kslope+b where m and b come from fit + ap_ufixed<18, 0> m = 0.03602; + ap_ufixed<18, 0> b = 0.66840; + int etaINT = int(round(double(m * track.kSlope() + b))); + int etaINTDisp = int(round(double(m * track.kSlopeAtMuon() + b))); double lsbEta = M_PI / pow(2, 12); int charge = 1; @@ -908,21 +1264,24 @@ void KMTFCore::setFourVectors(l1t::KMTFTrack& track) { //Also set PT =0 and dxy=0 if (track.stubs().size() == 1) { ptC = 0; - track.setCoordinatesAtMuon(track.curvatureAtMuon(), track.stubs()[0]->coord1(), track.phiBAtMuon()); - track.setCoordinatesAtVertex(track.curvatureAtVertex(), track.phiAtVertex(), 0); + track.setCoordinatesAtMuon( + track.curvatureAtMuon(), track.stubs()[0]->coord1(), track.phiBAtMuon(), track.zPosition(), track.kSlope()); + track.setCoordinatesAtVertex(track.curvatureAtVertex(), track.phiAtVertex(), 0, track.zPosition(), track.kSlope()); } track.setPt(ptC, ptU); //shift the dxy by 10 bits - track.setCoordinatesAtVertex(track.curvatureAtVertex(), track.phiAtVertex(), track.dxy() / 1024); + track.setCoordinatesAtVertex( + track.curvatureAtVertex(), track.phiAtVertex(), track.dxy() / 1024, track.zPosition(), track.kSlope()); //vertex double pt = ptC * 0.03125; double phi = (track.phiAtMuon() / 32) * M_PI / pow(2, 12); double eta = etaINT * lsbEta; + double etaDisp = etaINTDisp * lsbEta; track.setPtEtaPhi(pt, eta, phi); track.setCharge(charge); pt = double(ptLUT(track.curvatureAtMuon())) * 0.03125; - track.setPtEtaPhiDisplaced(pt, eta, phi); + track.setPtEtaPhiDisplaced(pt, etaDisp, phi); } bool KMTFCore::estimateChiSquare(l1t::KMTFTrack& track, bool vertex) { @@ -947,6 +1306,15 @@ bool KMTFCore::estimateChiSquare(l1t::KMTFTrack& track, bool vertex) { std::vector cut; const l1t::MuonStubRef& innerStub = track.stubs()[track.stubs().size() - 1]; + if (verbose_) { + edm::LogInfo("KMTFCore") << "[KMTF_CHI2DBG] begin " << (vertex ? "prompt" : "displaced") + << " chi2 fit: hitPattern=" << track.hitPattern() + << " thetaPattern=" << track.thetaDigiPattern() << " nStubs=" << track.stubs().size() + << " innerDepth=" << innerStub->depthRegion() << " curvMuon=" << track.curvatureAtMuon() + << " curvVtx=" << track.curvatureAtVertex() << " phiMuon=" << track.phiAtMuon() + << " phiVtx=" << track.phiAtVertex() << " z=" << track.zPosition() + << " kSlope=" << track.kSlope(); + } if (vertex) { K = track.curvatureAtVertex(); @@ -1019,14 +1387,19 @@ bool KMTFCore::estimateChiSquare(l1t::KMTFTrack& track, bool vertex) { chi = chi + absDelta; chiErr = chiErr + err; if (verbose_) { - edm::LogInfo("KMTFCore") << "Chi Square stub for track with pattern=" << track.hitPattern() - << " inner stub depth=" << innerStub->depthRegion() << "-> AK=" << int(AK) - << " stubDepth=" << stub->depthRegion() << " diff1=" << diffPhi << " diff2=" << diffPhiB - << " delta=" << absDelta << " absK=" << uint(absK) << " err=" << err; + edm::LogInfo("KMTFCore") << "[KMTF_CHI2DBG] per-stub " << (vertex ? "prompt" : "displaced") + << " innerDepth=" << innerStub->depthRegion() << " stubDepth=" << stub->depthRegion() + << " diffPhi=" << diffPhi << " diffPhiB=" << diffPhiB << " propC=" << float(propC) + << " Kshifted=" << int(Kshifted) << " AK=" << int(AK) << " deltaSigned=" << delta + << " deltaAbs=" << absDelta << " posErrA=" << positionError << " errBTerm=" + << uint(ap_fixed<8, 2>(propErrB[stub->depthRegion() - innerStub->depthRegion() - 1]) * + absK) + << " errTotal=" << err << " chiRunning=" << chi << " chiErrRunning=" << chiErr; } } if (verbose_) { - edm::LogInfo("KMTFCore") << "Chi Square =" << chi << " ChiSquare Error = " << chiErr; + edm::LogInfo("KMTFCore") << "[KMTF_CHI2DBG] final " << (vertex ? "prompt" : "displaced") << " chi=" << chi + << " chiErr=" << chiErr << " pass=" << (chi <= chiErr); } track.setApproxChi2(chi, chiErr, vertex); @@ -1073,6 +1446,12 @@ int KMTFCore::wrapAround(int value, int maximum) { return value; } +int KMTFCore::satKSlope(int k) { + const int kmax = (1 << (BITSKSLOPE - 1)) - 1; + const int kmin = -(1 << (BITSKSLOPE - 1)); + return k > kmax ? kmax : (k < kmin ? kmin : k); +} + int KMTFCore::encode(bool ownwheel, int sector, int tag) { int wheel = ownwheel ? 1 : 0; int phi = 0; @@ -1086,7 +1465,6 @@ int KMTFCore::encode(bool ownwheel, int sector, int tag) { std::pair KMTFCore::getByCode(const std::vector& tracks, int mask) { for (uint i = 0; i < tracks.size(); ++i) { - edm::LogInfo("KMTFCore") << "Code=" << tracks[i].hitPattern() << ", track=" << mask; if (tracks[i].hitPattern() == mask) return std::make_pair(true, i); } diff --git a/L1Trigger/Phase2L1GMT/src/L1TPhase2GMTBarrelStubProcessor.cc b/L1Trigger/Phase2L1GMT/src/L1TPhase2GMTBarrelStubProcessor.cc index b8743c7e58f5e..6da8b0b2b9f8b 100644 --- a/L1Trigger/Phase2L1GMT/src/L1TPhase2GMTBarrelStubProcessor.cc +++ b/L1Trigger/Phase2L1GMT/src/L1TPhase2GMTBarrelStubProcessor.cc @@ -114,53 +114,106 @@ l1t::MuonStub L1TPhase2GMTBarrelStubProcessor::buildStubNoEta(const L1Phase2MuDT return stub; } -l1t::MuonStubCollection L1TPhase2GMTBarrelStubProcessor::makeStubs(const L1Phase2MuDTPhContainer* phiContainer, - const L1MuDTChambThContainer* etaContainer) { +l1t::MuonStub L1TPhase2GMTBarrelStubProcessor::buildStubwithZandkSlope(const L1Phase2MuDTExtPhiThetaPair& pairs) { + const auto& phiS = pairs.phiDigi(); + static constexpr int ZCenterDigitizedPositive[] = {11730, 23327}; + static constexpr int ZCenterDigitizedNegative[] = {-11697, -23290}; + static constexpr int ZCenterDigitizedZero[] = {22}; + static constexpr float RadiusStationPhys[] = {445., 526., 635., 730.}; + int wheel = phiS.whNum(); + int abswheel = fabs(phiS.whNum()); + int sector = phiS.scNum(); + int station = phiS.stNum(); + + ap_uint<18> normalization0 = sector * ap_uint<15>(21845); + ap_int<18> normalization1 = ap_int<18>(ap_int<17>(phiS.phi()) * ap_ufixed<8, 0>(0.3183)); + ap_int<18> kmtf_phi = ap_int<18>(normalization0 + normalization1); + int phi = int(kmtf_phi); + float globalPhi = phi * M_PI / (1 << 17); + + // double globalPhi = (sector * 30) + phiS.phi() * 30. / 65535.; + int tag = phiS.index(); + + int bx = phiS.bxNum() - 20; + int quality = phiS.quality(); + uint tfLayer = phiS.stNum() - 1; + + // instantiate stub with eta1=0, eta2=0, etaQuality=0 values which eventually get written over with stub.setEta below. + l1t::MuonStub stub(wheel, sector, station, tfLayer, phi, phiS.phiBend(), tag, bx, quality, 0, 0, 0, 1); + + //defining z, k, zPhys, kPhys for case where theta digi exists + ap_int<16> z = pairs.thetaDigi().z(); + ap_int<16> k = pairs.thetaDigi().k(); + float zPhys = z * (1500. / (1 << 16)); + float kPhys = k * (2. / (1 << 16)); + + //slight asymmetry in wheels w.r.t. the origin calls for different z_center values for case where no theta digi was matched + //defining z_center adn k_center for when theta digi does NOT exist!! + int z_centerDigi; + if (wheel > 0) { + z_centerDigi = ZCenterDigitizedPositive[abswheel - 1]; + } else if (wheel < 0) { + z_centerDigi = ZCenterDigitizedNegative[abswheel - 1]; + } else { + z_centerDigi = ZCenterDigitizedZero[0]; + } + float z_centerPhys = z_centerDigi * (1500. / (1 << 16)); + float R_centerPhys = RadiusStationPhys[station - 1]; + float k_centerPhys = z_centerPhys / R_centerPhys; + int k_centerDigi = k_centerPhys * ((1 << 16) / 2.); + //very important + //clamp to a signed 16-bit range. upstream z and kSlope were declared with keyword int from phi-theta digi matched pair dataformat + //. this triggers only in wheels 1/2 wher k can exceed +/-1.0 for case of no theta digi!!!!!!!!!!!!! + //point back as close to the origin as we can. fixed range of kSlope binds the value + if (k_centerDigi > (1 << 15) - 1) { + k_centerDigi = (1 << 15) - 1; + } + if (k_centerDigi < -(1 << 15)) { + k_centerDigi = -(1 << 15); + } + // check if theta digi exists with non-default constructor quality --> use z, k with etaQuality=3 from theta digi. + // if theta digi has no real data, use z_center and slope which points to origin with etaQuality=0.. + // stub set to etaQuality==3 if theta digi exists, 0 if not. + if (pairs.thetaDigi().quality() >= 0) { + stub.setEta(z, k, 3); + stub.setOfflineQuantities(globalPhi, float(phiS.phiBend() * 0.49e-3), zPhys, kPhys); + } else { + stub.setEta(z_centerDigi, k_centerDigi, 0); + stub.setOfflineQuantities(globalPhi, float(phiS.phiBend() * 0.49e-3), z_centerPhys, k_centerPhys); + } + return stub; +} + +l1t::MuonStubCollection L1TPhase2GMTBarrelStubProcessor::makeStubs( + const L1Phase2MuDTExtPhiThetaPairContainer* pairContainer) { l1t::MuonStubCollection out; for (int bx = minBX_; bx <= maxBX_; bx++) { ostringstream os; if (verbose_ == 2) os << "PATTERN "; - for (int wheel = -2; wheel <= 2; wheel++) { - for (int sector = 0; sector < 12; sector++) { - for (int station = 1; station < 5; station++) { - bool hasEta = false; - const L1MuDTChambThDigi* tseta = etaContainer->chThetaSegm(wheel, station, sector, bx); - if (tseta != nullptr) { - hasEta = true; - } - - for (const auto& phiDigi : *phiContainer->getContainer()) { - if ((phiDigi.bxNum() - 20) != bx || phiDigi.whNum() != wheel || phiDigi.scNum() != sector || - phiDigi.stNum() != station) - continue; - if (phiDigi.quality() < minPhiQuality_) - continue; - - if (verbose_ == 2) { - ap_uint<64> wphi = ap_uint<17>(phiDigi.phi()); - ap_uint<64> wphib = ap_uint<13>(phiDigi.phiBend()); - ap_uint<64> wr1 = ap_uint<21>(0); - ap_uint<64> wq = ap_uint<4>(phiDigi.quality()); - ap_uint<64> wr2 = ap_uint<9>(0); - ap_uint<64> sN = 0; - sN = sN | wphi; - sN = sN | (wphib << 17); - sN = sN | (wr1 << 30); - sN = sN | (wq << 51); - sN = sN | (wr2 << 55); - os << std::setw(0) << std::dec << sector << " " << wheel << " " << station << " "; - os << std::uppercase << std::setfill('0') << std::setw(16) << std::hex << uint64_t(sN) << " "; - } - - if (hasEta) { - out.push_back(buildStub(phiDigi, tseta)); - } else { - out.push_back(buildStubNoEta(phiDigi)); - } - } - } + for (const auto& pair : pairContainer->getContainer()) { + const auto& phiDigi = pair.phiDigi(); + if ((phiDigi.bxNum() - 20) != bx) + continue; + if (phiDigi.quality() < minPhiQuality_) + continue; + + if (verbose_ == 2) { + ap_uint<64> wphi = ap_uint<17>(phiDigi.phi()); + ap_uint<64> wphib = ap_uint<13>(phiDigi.phiBend()); + ap_uint<64> wr1 = ap_uint<21>(0); + ap_uint<64> wq = ap_uint<4>(phiDigi.quality()); + ap_uint<64> wr2 = ap_uint<9>(0); + ap_uint<64> sN = 0; + sN = sN | wphi; + sN = sN | (wphib << 17); + sN = sN | (wr1 << 30); + sN = sN | (wq << 51); + sN = sN | (wr2 << 55); + os << std::setw(0) << std::dec << phiDigi.scNum() << " " << phiDigi.whNum() << " " << phiDigi.stNum() << " "; + os << std::uppercase << std::setfill('0') << std::setw(16) << std::hex << uint64_t(sN) << " "; } + out.push_back(buildStubwithZandkSlope(pair)); } if (verbose_ == 2) edm::LogInfo("BarrelStub") << os.str() << std::endl;