Partially fix CLCT position bias after CCLUT (CCLUT-8)

DataFormats/CSCDigi/interface/CSCCLCTDigi.h

@@ -74,13 +74,18 @@ class CSCCLCTDigi {
   /// set the slope
   void setSlope(const uint16_t slope);
 
+  /// slope in number of half-strips/layer
+  float getFractionalSlope(const uint16_t slope = 5) const;
+
   /// return striptype
   uint16_t getStripType() const { return striptype_; }
 
   /// set stripType
   void setStripType(const uint16_t stripType) { striptype_ = stripType; }
 
-  /// return bend (left or right)
+  /// return bending
+  /// 0: left-bending (negative delta-strip)
+  /// 1: right-bending (positive delta-strip)
   uint16_t getBend() const { return bend_; }
 
   /// set bend

DataFormats/CSCDigi/interface/CSCCorrelatedLCTDigi.h

@@ -21,7 +21,7 @@ class CSCCorrelatedLCTDigi {
   enum LCTKeyStripMasks { kEightStripMask = 0x1, kQuartStripMask = 0x1, kHalfStripMask = 0xff };
   enum LCTKeyStripShifts { kEightStripShift = 9, kQuartStripShift = 8, kHalfStripShift = 0 };
   // temporary to facilitate CCLUT-EMTF/OMTF integration studies
-  enum LCTPatternMasks { kRun3SlopeMask = 0x1f, kRun3PatternMask = 0x7, kLegacyPatternMask = 0xf };
+  enum LCTPatternMasks { kRun3SlopeMask = 0xf, kRun3PatternMask = 0x7, kLegacyPatternMask = 0xf };
   enum LCTPatternShifts { kRun3SlopeShift = 7, kRun3PatternShift = 4, kLegacyPatternShift = 0 };
   enum class Version { Legacy = 0, Run3 };
 

DataFormats/CSCDigi/src/CSCCLCTDigi.cc

@@ -115,6 +115,21 @@ void CSCCLCTDigi::setSlope(const uint16_t slope) {
   setDataWord(slope, pattern_, kRun3SlopeShift, kRun3SlopeMask);
 }
 
+// slope in number of half-strips/layer
+float CSCCLCTDigi::getFractionalSlope(const uint16_t nBits) const {
+  if (isRun3()) {
+    const float minSlope = 0;
+    const float maxSlope = 2.5;
+    const int range = pow(2, nBits);
+    const float deltaSlope = (maxSlope - minSlope) / range;
+    const float slopeValue = minSlope + deltaSlope * getSlope();
+    return (2 * getBend() - 1) * slopeValue;
+  } else {
+    int slope[11] = {0, 0, -8, 8, -6, 6, -4, 4, -2, 2, 0};
+    return float(slope[getPattern()] / 5.);
+  }
+}
+
 uint16_t CSCCLCTDigi::getKeyStrip(const uint16_t n) const {
   // 10-bit case for strip data word
   if (compCode_ != -1 and n == 8) {
@@ -132,9 +147,9 @@ uint16_t CSCCLCTDigi::getKeyStrip(const uint16_t n) const {
 
 /// return the fractional strip (middle of the strip)
 float CSCCLCTDigi::getFractionalStrip(const uint16_t n) const {
-  if (n == 8) {
+  if (compCode_ != -1 and n == 8) {
     return 0.125f * (getKeyStrip(n) + 0.5);
-  } else if (n == 4) {
+  } else if (compCode_ != -1 and n == 4) {
     return 0.25f * (getKeyStrip(n) + 0.5);
   } else {
     return 0.5f * (getKeyStrip(n) + 0.5);

L1Trigger/CSCTriggerPrimitives/interface/CSCCathodeLCTProcessor.h

@@ -158,7 +158,7 @@ class CSCCathodeLCTProcessor : public CSCBaseboard {
   int calculateComparatorCode(const std::array<std::array<int, 3>, 6>& halfStripPattern) const;
 
   // sets the 1/4 and 1/8 strip bits given a floating point position offset
-  void assignPositionCC(const unsigned offset, uint16_t& halfstrip, bool& quartstrip, bool& eightstrip) const;
+  void assignPositionCC(const unsigned offset, std::tuple<uint16_t, bool, bool>& returnValue) const;
 
   // runs the CCLUT procedure
   void runCCLUT(CSCCLCTDigi& digi) const;

L1Trigger/CSCTriggerPrimitives/src/CSCCathodeLCTProcessor.cc

@@ -3,7 +3,6 @@
 #include <iomanip>
 #include <iostream>
 #include <memory>
-#include <bitset>
 
 // Default values of configuration parameters.
 const unsigned int CSCCathodeLCTProcessor::def_fifo_tbins = 12;
@@ -1310,9 +1309,7 @@ int CSCCathodeLCTProcessor::calculateComparatorCode(const std::array<std::array<
 }
 
 void CSCCathodeLCTProcessor::assignPositionCC(const unsigned offset,
-                                              uint16_t& halfstrip,
-                                              bool& quartstrip,
-                                              bool& eightstrip) const {
+                                              std::tuple<uint16_t, bool, bool>& returnValue) const {
   /*
     | Value | Half-Strip Offset  | Delta Half-Strip  | Quarter-Strip Bit  | Eighth-Strip Bit |
     |-------|--------------------|-------------------|--------------------|------------------|
@@ -1333,17 +1330,25 @@ void CSCCathodeLCTProcessor::assignPositionCC(const unsigned offset,
     |   14  |   7/4              |   1               |   1                |   1              |
     |   15  |   2                |   2               |   0                |   0              |
    */
-  if (offset <= 2)
-    halfstrip -= 2;
-  else if (offset > 2 and offset <= 6)
-    halfstrip--;
-  else if (offset > 10 and offset <= 14)
-    halfstrip += 1;
-  else if (offset == 15)
-    halfstrip += 2;
-
-  quartstrip = std::bitset<4>(offset + 1)[1];
-  eightstrip = !std::bitset<4>(offset)[0];
+  std::vector<std::tuple<uint16_t, bool, bool>> my_tuple = {
+      {-2, false, true},
+      {-2, true, false},
+      {-2, true, true},
+      {-1, false, false},
+      {-1, false, true},
+      {-1, true, false},
+      {-1, true, true},
+      {0, false, false},
+      {0, false, true},
+      {0, true, false},
+      {0, true, true},
+      {1, false, false},
+      {1, false, true},
+      {1, true, false},
+      {1, true, true},
+      {2, false, false},
+  };
+  returnValue = my_tuple[offset];
 }
 
 void CSCCathodeLCTProcessor::runCCLUT(CSCCLCTDigi& digi) const {
@@ -1403,24 +1408,23 @@ void CSCCathodeLCTProcessor::runCCLUT(CSCCLCTDigi& digi) const {
   unsigned run2PatternCC(lutpatconv_[pattern]->lookup(comparatorCode));
 
   // if the slope is negative, set bending to 0
-  if (slopeCC < 16)
-    digi.setBend(0);
-  else
-    digi.setBend(1);
+  const bool slopeCCSign((slopeCC >> 4) & 0x1);
+  const unsigned slopeCCValue(slopeCC & 0xf);
+  digi.setBend(slopeCCSign);
 
   // calculate the new position
   uint16_t halfstrip = digi.getKeyStrip();
-  bool quartstrip = false;
-  bool eightstrip = false;
-  assignPositionCC(positionCC, halfstrip, quartstrip, eightstrip);
+  std::tuple<uint16_t, bool, bool> halfstripoffset;
+  assignPositionCC(positionCC, halfstripoffset);
+  halfstrip += std::get<0>(halfstripoffset);
 
   // store the new 1/2, 1/4 and 1/8 strip positions
   digi.setStrip(halfstrip - digi.getCFEB() * 32);
-  digi.setQuartStrip(quartstrip);
-  digi.setEightStrip(eightstrip);
+  digi.setQuartStrip(std::get<1>(halfstripoffset));
+  digi.setEightStrip(std::get<2>(halfstripoffset));
 
   // store the bending angle value in the pattern data member
-  digi.setSlope(slopeCC);
+  digi.setSlope(slopeCCValue);
 
   // set the quasi Run-2 pattern - to accommodate integration with EMTF/OMTF
   if (run2PatternCC == 0) {
@@ -1444,7 +1448,6 @@ void CSCCathodeLCTProcessor::runCCLUT(CSCCLCTDigi& digi) const {
 }
 
 unsigned CSCCathodeLCTProcessor::convertSlopeToRun2Pattern(const unsigned slope) const {
-  // interesting: the "right-bending" patterns have a negative slope
   const unsigned slopeList[32] = {2,  2,  2, 4, 4, 4, 6, 6, 6, 6, 8, 8, 8, 8, 10, 10,
                                   10, 10, 9, 9, 9, 9, 7, 7, 7, 7, 5, 5, 5, 3, 3,  3};
   return slopeList[slope];

L1Trigger/CSCTriggerPrimitives/src/CSCMotherboard.cc

@@ -489,9 +489,6 @@ CSCCorrelatedLCTDigi CSCMotherboard::constructLCTs(const CSCALCTDigi& aLCT,
                                                    int trknmb) const {
   // CLCT pattern number
   unsigned int pattern = encodePattern(cLCT.getPattern());
-  if (use_run3_patterns_ and use_comparator_codes_) {
-    pattern = cLCT.getSlope();
-  }
 
   // LCT quality number
   unsigned int quality;

L1Trigger/CSCTriggerPrimitives/test/macros/CCLUTLinearFitWriter.cpp

@@ -130,7 +130,8 @@ void writeHeaderPosOffsetLUT(ofstream& file);
 void writeHeaderSlopeLUT(ofstream& file);
 unsigned firmwareWord(const unsigned quality, const unsigned slope, const unsigned offset);
 void setDataWord(unsigned& word, const unsigned newWord, const unsigned shift, const unsigned mask);
-unsigned assign(const float fvalue, const float fmin, const float fmax, const unsigned nbits);
+unsigned assignPosition(const float fvalue, const float fmin, const float fmax, const unsigned nbits);
+unsigned assignBending(const float fvalue, const float fmin, const float fmax, const unsigned nbits);
 
 int CCLUTLinearFitWriter(unsigned N_LAYER_REQUIREMENT = 3) {
   //all the patterns we will fit
@@ -340,10 +341,16 @@ int CCLUTLinearFitWriter(unsigned N_LAYER_REQUIREMENT = 3) {
       const float fmaxOffset = 2;
       const float fminOffset = -1.75;
       const float fmaxSlope = 2.5;
-      const float fminSlope = -2.5;
+      const float fminSlope = 0;
 
-      const unsigned offset_bin = assign(offset, fminOffset, fmaxOffset, 4);
-      const unsigned slope_bin = assign(slope, fminSlope, fmaxSlope, 5);
+      // negative bending -> 0
+      // positive bending -> 1
+      const bool slope_sign(slope >= 0);
+
+      const unsigned offset_bin = assignPosition(offset, fminOffset, fmaxOffset, 4);
+      unsigned slope_bin = assignBending(std::abs(slope), fminSlope, fmaxSlope, 4);
+      if (slope_sign)
+        slope_bin += 16;
       const unsigned fwword = firmwareWord(0, slope_bin, offset_bin);
 
       // write to output files
@@ -496,7 +503,28 @@ void writeHeaderSlopeLUT(ofstream& file) {
        << "#<header> v1.0 12 32 </header>\n";
 }
 
-unsigned assign(const float fvalue, const float fmin, const float fmax, const unsigned nbits) {
+unsigned assignPosition(const float fvalue, const float fmin, const float fmax, const unsigned nbits) {
+  bool debug;
+  unsigned value = 0;
+  const unsigned range = pow(2, nbits);
+  const unsigned minValue = 0;
+  const unsigned maxValue = range - 1;
+  const double fdelta = (fmax - fmin) / range;
+
+  if (fvalue >= fmax) {
+    value = maxValue;
+  } else if (fvalue <= fmin) {
+    value = minValue;
+  } else {
+    value = std::min(unsigned(std::ceil((fvalue - fmin) / fdelta)), maxValue);
+  }
+  if (debug)
+    std::cout << "fvalue " << fvalue << " " << fmin << " " << fmax << " " << nbits << " " << value << std::endl;
+
+  return value;
+}
+
+unsigned assignBending(const float fvalue, const float fmin, const float fmax, const unsigned nbits) {
   bool debug;
   unsigned value = 0;
   const unsigned range = pow(2, nbits);
@@ -509,7 +537,7 @@ unsigned assign(const float fvalue, const float fmin, const float fmax, const un
   } else if (fvalue <= fmin) {
     value = minValue;
   } else {
-    value = int(std::floor((fvalue - fmin) / fdelta));
+    value = std::min(unsigned(std::floor((fvalue - fmin) / fdelta)), maxValue);
   }
   if (debug)
     std::cout << "fvalue " << fvalue << " " << fmin << " " << fmax << " " << nbits << " " << value << std::endl;
@@ -520,11 +548,12 @@ unsigned assign(const float fvalue, const float fmin, const float fmax, const un
 unsigned firmwareWord(const unsigned quality, const unsigned slope, const unsigned offset) {
   /* construct fw dataword:
      [8:0] is quality (set all to 0 for now)
-     [13, 9] is slope
+     [12:9] is slope value
+     [13] is slope sign
      [17:14] is offset
   */
   enum Masks { OffsetMask = 0xf, SlopeMask = 0x1f, QualityMask = 0x1ff };
-  enum Shifts { OffsetShift = 13, SlopeShift = 9, QualityShift = 0 };
+  enum Shifts { OffsetShift = 14, SlopeShift = 9, QualityShift = 0 };
 
   unsigned fwword = 0;
   setDataWord(fwword, quality, QualityShift, QualityMask);

Validation/MuonGEMDigis/src/GEMDigiMatcher.cc

@@ -13,7 +13,7 @@ GEMDigiMatcher::GEMDigiMatcher(const edm::ParameterSet& pset, edm::ConsumesColle
   maxBXDigi_ = gemDigi.getParameter<int>("maxBX");
   matchDeltaStrip_ = gemDigi.getParameter<int>("matchDeltaStrip");
   verboseDigi_ = gemDigi.getParameter<int>("verbose");
-  matchToSimLink_ = gemDigi.getParameter<int>("matchToSimLink");
+  matchToSimLink_ = gemDigi.getParameter<bool>("matchToSimLink");
 
   const auto& gemPad = pset.getParameterSet("gemPadDigi");
   minBXPad_ = gemPad.getParameter<int>("minBX");
@@ -33,7 +33,9 @@ GEMDigiMatcher::GEMDigiMatcher(const edm::ParameterSet& pset, edm::ConsumesColle
   // make a new simhits matcher
   muonSimHitMatcher_.reset(new GEMSimHitMatcher(pset, std::move(iC)));
 
-  gemSimLinkToken_ = iC.consumes<edm::DetSetVector<GEMDigiSimLink>>(gemSimLink.getParameter<edm::InputTag>("inputTag"));
+  if (matchToSimLink_)
+    gemSimLinkToken_ =
+        iC.consumes<edm::DetSetVector<GEMDigiSimLink>>(gemSimLink.getParameter<edm::InputTag>("inputTag"));
   gemDigiToken_ = iC.consumes<GEMDigiCollection>(gemDigi.getParameter<edm::InputTag>("inputTag"));
   gemPadToken_ = iC.consumes<GEMPadDigiCollection>(gemPad.getParameter<edm::InputTag>("inputTag"));
   gemClusterToken_ = iC.consumes<GEMPadDigiClusterCollection>(gemCluster.getParameter<edm::InputTag>("inputTag"));
@@ -45,7 +47,8 @@ GEMDigiMatcher::GEMDigiMatcher(const edm::ParameterSet& pset, edm::ConsumesColle
 void GEMDigiMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
   muonSimHitMatcher_->init(iEvent, iSetup);
 
-  iEvent.getByToken(gemSimLinkToken_, gemDigisSLH_);
+  if (matchToSimLink_)
+    iEvent.getByToken(gemSimLinkToken_, gemDigisSLH_);
   iEvent.getByToken(gemDigiToken_, gemDigisH_);
   iEvent.getByToken(gemPadToken_, gemPadsH_);
   iEvent.getByToken(gemClusterToken_, gemClustersH_);
@@ -73,7 +76,8 @@ void GEMDigiMatcher::match(const SimTrack& t, const SimVertex& v) {
     return;
 
   // now match the digis
-  matchDigisSLToSimTrack(gemDigisSL);
+  if (matchToSimLink_)
+    matchDigisSLToSimTrack(gemDigisSL);
   matchDigisToSimTrack(gemDigis);
   matchPadsToSimTrack(gemPads);
   matchClustersToSimTrack(gemClusters);

Validation/MuonHits/interface/CSCSimHitMatcher.h

@@ -55,6 +55,7 @@ class CSCSimHitMatcher : public MuonSimHitMatcher {
 
   // local bending in a CSC chamber
   float LocalBendingInChamber(unsigned int detid) const;
+  void fitHitsInChamber(unsigned int detid, float& mean, float& slope) const;
 
   // calculate average strip number for a provided collection of simhits
   float simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const;

Validation/MuonHits/interface/MuonSimHitMatcher.h

@@ -66,6 +66,8 @@ class MuonSimHitMatcher {
   // calculate the average position at the second station
   GlobalPoint simHitsMeanPositionStation(int n) const;
 
+  const TrackingGeometry* geometry() { return geometry_; }
+
 protected:
   std::vector<unsigned int> getIdsOfSimTrackShower(unsigned trk_id,
                                                    const edm::SimTrackContainer& simTracks,

Validation/MuonHits/src/CSCSimHitMatcher.cc

@@ -1,4 +1,6 @@
 #include "Validation/MuonHits/interface/CSCSimHitMatcher.h"
+#include "TGraphErrors.h"
+#include "TF1.h"
 
 using namespace std;
 
@@ -204,6 +206,43 @@ float CSCSimHitMatcher::LocalBendingInChamber(unsigned int detid) const {
   return deltaPhi(phi_layer6, phi_layer1);
 }
 
+// difference in strip per layer
+void CSCSimHitMatcher::fitHitsInChamber(unsigned int detid, float& intercept, float& slope) const {
+  const CSCDetId cscid(detid);
+
+  const auto& sim_hits = hitsInChamber(detid);
+
+  if (sim_hits.empty())
+    return;
+
+  vector<float> x;
+  vector<float> y;
+  vector<float> xe;
+  vector<float> ye;
+
+  const float HALF_STRIP_ERROR = 0.288675;
+
+  for (const auto& h : sim_hits) {
+    const LocalPoint& lp = h.entryPoint();
+    const auto& d = h.detUnitId();
+    float s = dynamic_cast<const CSCGeometry*>(geometry_)->layer(d)->geometry()->strip(lp);
+    // shift to key half strip layer (layer 3)
+    x.push_back(CSCDetId(d).layer() - 3);
+    y.push_back(s);
+    xe.push_back(float(0));
+    ye.push_back(2 * HALF_STRIP_ERROR);
+  }
+  if (x.size() < 2)
+    return;
+
+  std::unique_ptr<TGraphErrors> gr(new TGraphErrors(x.size(), &x[0], &y[0], &xe[0], &ye[0]));
+  std::unique_ptr<TF1> fit(new TF1("fit", "pol1", -3, 4));
+  gr->Fit("fit", "EMQ");
+
+  intercept = fit->GetParameter(0);
+  slope = fit->GetParameter(1);
+}
+
 float CSCSimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
   if (sim_hits.empty())
     return -1.f;