The effect of module sum split over multiple towers on physics objects

DataFormats/L1THGCal/interface/HGCalTowerID.h

@@ -26,15 +26,16 @@ namespace l1t {
 
     unsigned short rawId() const { return rawId_; }
 
-  private:
-    uint32_t rawId_;
     static const int subDetMask = 0x1;  // two for now 0 is HGC and 1 is HFNose
     static const int subDetShift = 16;
     static const int zsideMask = 0x1;
     static const int zsideShift = 15;
     static const int coordMask = 0x007F;
     static const int coord1Shift = 7;
     static const int coord2Shift = 0;
+
+  private:
+    uint32_t rawId_;
   };
 
   struct HGCalTowerCoord {

DataFormats/L1THGCal/interface/HGCalTowerMap.h

@@ -22,7 +22,7 @@ namespace l1t {
 
     const HGCalTowerMap& operator+=(const HGCalTowerMap& map);
 
-    bool addEt(short bin_id, float etEm, float etHad);
+    bool addEt(const std::unordered_map<unsigned short, float>& towerIDandShares, float etEm, float etHad);
 
     unsigned nTowers() const { return towerMap_.size(); }
     const std::unordered_map<unsigned short, l1t::HGCalTower>& towers() const { return towerMap_; }

DataFormats/L1THGCal/src/HGCalTowerMap.cc

@@ -30,12 +30,13 @@ const HGCalTowerMap& HGCalTowerMap::operator+=(const HGCalTowerMap& map) {
   return *this;
 }
 
-bool HGCalTowerMap::addEt(short bin_id, float etEm, float etHad) {
-  auto this_tower = towerMap_.find(bin_id);
-  if (this_tower == towerMap_.end())
-    return false;
-  this_tower->second.addEtEm(etEm);
-  this_tower->second.addEtHad(etHad);
-
+bool HGCalTowerMap::addEt(const std::unordered_map<unsigned short, float>& towerIDandShares, float etEm, float etHad) {
+  for (const auto& towerIDandShare : towerIDandShares) {
+    auto this_tower = towerMap_.find(towerIDandShare.first);
+    if (this_tower == towerMap_.end())
+      return false;
+    this_tower->second.addEtEm(etEm * towerIDandShare.second);
+    this_tower->second.addEtHad(etHad * towerIDandShare.second);
+  }
   return true;
 }

L1Trigger/L1THGCal/interface/HGCalTriggerTowerGeometryHelper.h

@@ -15,6 +15,7 @@
 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
 #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"
 #include "DataFormats/L1THGCal/interface/HGCalTriggerSums.h"
+#include "DataFormats/ForwardDetId/interface/HGCalTriggerModuleDetId.h"
 
 #include <vector>
 #include <unordered_map>
@@ -32,17 +33,29 @@ class HGCalTriggerTowerGeometryHelper {
 
   void eventSetup(const edm::EventSetup& es) { triggerTools_.eventSetup(es); }
 
+  unsigned packLayerSubdetWaferId(int subdet, int layer, int moduleU, int moduleV) const;
+  unsigned packTowerIDandShare(int towerEta, int towerPhi, int towerShare) const;
+  void unpackTowerIDandShare(unsigned towerIDandShare, int& towerEta_raw, int& towerPhi_raw, int& towerShare) const;
+  int moveToCorrectSector(int towerPhi_raw, int sector) const;
+  void reverseXaxis(int& towerPhi) const;
+
   const std::vector<l1t::HGCalTowerCoord>& getTowerCoordinates() const;
 
   unsigned short getTriggerTowerFromEtaPhi(const float& eta, const float& phi) const;
-  unsigned short getTriggerTower(const l1t::HGCalTriggerCell&) const;
-  unsigned short getTriggerTower(const l1t::HGCalTriggerSums&) const;
+  std::unordered_map<unsigned short, float> getTriggerTower(const l1t::HGCalTriggerCell&) const;
+  std::unordered_map<unsigned short, float> getTriggerTower(const l1t::HGCalTriggerSums&) const;
 
   const bool isNose() { return doNose_; }
 
 private:
+  static const int towerShareMask = 0x7F;
+  static const int towerShareShift = 14;
+  static const int signMask = 0x1;
+  static const int sign1Shift = 21;
+  static const int sign2Shift = 22;
   std::vector<l1t::HGCalTowerCoord> tower_coords_;
   std::unordered_map<unsigned, short> cells_to_trigger_towers_;
+  std::unordered_map<unsigned, std::vector<unsigned>> modules_to_trigger_towers_;
 
   bool doNose_;
   double minEta_;
@@ -55,6 +68,13 @@ class HGCalTriggerTowerGeometryHelper {
   std::vector<double> binsEta_;
   std::vector<double> binsPhi_;
 
+  bool splitModuleSum_;
+  int splitDivisorSilic_;
+  int splitDivisorScint_;
+  int rotate180Deg_;
+  int rotate120Deg_;
+  int reverseX_;
+
   HGCalTriggerTools triggerTools_;
 };
 

L1Trigger/L1THGCal/python/customTowers.py

@@ -1,4 +1,5 @@
 import FWCore.ParameterSet.Config as cms
+from L1Trigger.L1THGCal.hgcalTowerMapProducer_cfi import L1TTriggerTowerConfig_energySplit
 import math
 
 def custom_towers_unclustered_tc(process):
@@ -35,6 +36,10 @@ def custom_towers_etaphi(process,
     parameters_towers_2d.L1TTriggerTowerConfig.binsPhi = cms.vdouble(binsPhi)
     return process
 
+def custom_towers_energySplit(process):
+    parameters_towers_2d = L1TTriggerTowerConfig_energySplit.clone()
+    process.hgcalTowerMapProducer.ProcessorParameters.towermap_parameters.L1TTriggerTowerConfig = parameters_towers_2d
+    return process
 
 def custom_towers_map(process,
         towermapping='L1Trigger/L1THGCal/data/tower_mapping_hgcroc_eta-phi_v3.txt',

L1Trigger/L1THGCal/python/hgcalTowerMapProducer_cfi.py

@@ -10,7 +10,21 @@
                                         nBinsEta=cms.int32(18),
                                         nBinsPhi=cms.int32(72),
                                         binsEta=cms.vdouble(),
-                                        binsPhi=cms.vdouble())
+                                        binsPhi=cms.vdouble(),
+                                        splitModuleSum=cms.bool(False))
+
+L1TTriggerTowerConfig_energySplit = cms.PSet(readMappingFile=cms.bool(False),
+                                        doNose=cms.bool(False),
+                                        minEta=cms.double(1.305),
+                                        maxEta=cms.double(3.045),
+                                        minPhi=cms.double(-1*math.pi),
+                                        maxPhi=cms.double(math.pi),
+                                        nBinsEta=cms.int32(20),
+                                        nBinsPhi=cms.int32(72),
+                                        binsEta=cms.vdouble(),
+                                        binsPhi=cms.vdouble(),
+                                        splitModuleSum=cms.bool(True),
+                                        moduleTowerMapping=cms.FileInPath("L1Trigger/L1THGCal/data/tower_per_module_silic8_scint16.txt"))
 
 towerMap2D_parValues = cms.PSet( useLayerWeights = cms.bool(False),
                                  layerWeights = cms.vdouble(),

L1Trigger/L1THGCal/src/HGCalTriggerTowerGeometryHelper.cc

@@ -19,7 +19,8 @@ HGCalTriggerTowerGeometryHelper::HGCalTriggerTowerGeometryHelper(const edm::Para
       nBinsEta_(conf.getParameter<int>("nBinsEta")),
       nBinsPhi_(conf.getParameter<int>("nBinsPhi")),
       binsEta_(conf.getParameter<std::vector<double> >("binsEta")),
-      binsPhi_(conf.getParameter<std::vector<double> >("binsPhi")) {
+      binsPhi_(conf.getParameter<std::vector<double> >("binsPhi")),
+      splitModuleSum_(conf.getParameter<bool>("splitModuleSum")) {
   if (!binsEta_.empty() && ((unsigned int)(binsEta_.size()) != nBinsEta_ + 1)) {
     throw edm::Exception(edm::errors::Configuration, "Configuration")
         << "HGCalTriggerTowerGeometryHelper nBinsEta for the tower map not consistent with binsEta size" << std::endl;
@@ -79,6 +80,129 @@ HGCalTriggerTowerGeometryHelper::HGCalTriggerTowerGeometryHelper(const edm::Para
     }
     l1tTriggerTowerMappingStream.close();
   }
+
+  if (splitModuleSum_) {
+    //variables for transforming towers
+    rotate180Deg_ = int(nBinsPhi_) / 2;
+    rotate120Deg_ = int(nBinsPhi_) / 3;
+    reverseX_ = int(nBinsPhi_) / 2 - 1;
+
+    std::ifstream moduleTowerMappingStream(conf.getParameter<edm::FileInPath>("moduleTowerMapping").fullPath());
+    if (!moduleTowerMappingStream.is_open()) {
+      throw cms::Exception("MissingDataFile") << "Cannot open HGCalTowerMapProducer moduleTowerMapping file\n";
+    }
+    //get split divisors
+    std::string line;
+    std::getline(moduleTowerMappingStream, line);  //Skip row
+    std::getline(moduleTowerMappingStream, line);
+    std::stringstream ss(line);
+    ss >> splitDivisorSilic_ >> splitDivisorScint_;
+
+    //get towers and module sum shares
+    std::getline(moduleTowerMappingStream, line);  //Skip row
+    std::getline(moduleTowerMappingStream, line);  //Skip row
+    const int minNumOfWordsPerRow = 5;
+    const int numOfWordsPerTower = 3;
+    for (std::string line; std::getline(moduleTowerMappingStream, line);) {
+      int numOfWordsInThisRow = 0;
+      for (std::string::size_type i = 0; i < line.size(); i++) {
+        if (line[i] != ' ' && line[i + 1] == ' ') {
+          numOfWordsInThisRow++;
+        }
+      }
+      if (numOfWordsInThisRow < minNumOfWordsPerRow) {
+        throw edm::Exception(edm::errors::Configuration, "Configuration")
+            << "HGCalTriggerTowerGeometryHelper warning: Incorrect/incomplete values for module ID in the mapping "
+               "file.\n"
+            << "The incorrect line is:" << line << std::endl;
+      }
+      int subdet;
+      int layer;
+      int moduleU;
+      int moduleV;
+      int numTowers;
+      std::stringstream ss(line);
+      ss >> subdet >> layer >> moduleU >> moduleV >> numTowers;
+      if (numOfWordsInThisRow != (numTowers * numOfWordsPerTower + minNumOfWordsPerRow)) {
+        throw edm::Exception(edm::errors::Configuration, "Configuration")
+            << "HGCalTriggerTowerGeometryHelper warning: Incorrect/incomplete values for module ID or tower "
+               "share/eta/phi in the mapping file.\n"
+            << "The incorrect line is:" << line << std::endl;
+      }
+      unsigned packed_modID = packLayerSubdetWaferId(subdet, layer, moduleU, moduleV);
+      std::vector<unsigned> towers;
+      for (int itr_tower = 0; itr_tower < numTowers; itr_tower++) {
+        int iEta_raw;
+        int iPhi_raw;
+        int towerShare;
+        ss >> iEta_raw >> iPhi_raw >> towerShare;
+        int splitDivisor = (subdet == 2) ? splitDivisorScint_ : splitDivisorSilic_;
+        if ((towerShare > splitDivisor) || (towerShare < 1)) {
+          throw edm::Exception(edm::errors::Configuration, "Configuration")
+              << "HGCalTriggerTowerGeometryHelper warning: invalid tower share in the mapping file.\n"
+              << "Tower share must be a positive integer and less than splitDivisor. The incorrect values found for "
+                 "module ID:"
+              << std::endl
+              << "subdet=" << subdet << ", l=" << layer << ", u=" << moduleU << ", v=" << moduleV << std::endl;
+        }
+        towers.push_back(packTowerIDandShare(iEta_raw, iPhi_raw, towerShare));
+      }
+      modules_to_trigger_towers_[packed_modID] = towers;
+    }
+    moduleTowerMappingStream.close();
+  }
+}
+
+unsigned HGCalTriggerTowerGeometryHelper::packLayerSubdetWaferId(int subdet, int layer, int moduleU, int moduleV) const {
+  unsigned packed_modID = 0;
+  packed_modID |= ((subdet & HGCalTriggerModuleDetId::kHGCalTriggerSubdetMask)
+                   << HGCalTriggerModuleDetId::kHGCalTriggerSubdetOffset);
+  packed_modID |= ((layer & HGCalTriggerModuleDetId::kHGCalLayerMask) << HGCalTriggerModuleDetId::kHGCalLayerOffset);
+  packed_modID |=
+      ((moduleU & HGCalTriggerModuleDetId::kHGCalModuleUMask) << HGCalTriggerModuleDetId::kHGCalModuleUOffset);
+  packed_modID |=
+      ((moduleV & HGCalTriggerModuleDetId::kHGCalModuleVMask) << HGCalTriggerModuleDetId::kHGCalModuleVOffset);
+  return packed_modID;
+}
+
+unsigned HGCalTriggerTowerGeometryHelper::packTowerIDandShare(int iEta_raw, int iPhi_raw, int towerShare) const {
+  unsigned packed_towerIDandShare = 0;
+  unsigned iEtaAbs = std::abs(iEta_raw);
+  unsigned iEtaSign = std::signbit(iEta_raw);
+  unsigned iPhiAbs = std::abs(iPhi_raw);
+  unsigned iPhiSign = std::signbit(iPhi_raw);
+  packed_towerIDandShare |= ((iEtaAbs & l1t::HGCalTowerID::coordMask) << l1t::HGCalTowerID::coord1Shift);
+  packed_towerIDandShare |= ((iEtaSign & signMask) << sign1Shift);
+  packed_towerIDandShare |= ((iPhiAbs & l1t::HGCalTowerID::coordMask) << l1t::HGCalTowerID::coord2Shift);
+  packed_towerIDandShare |= ((iPhiSign & signMask) << sign2Shift);
+  packed_towerIDandShare |= ((towerShare & towerShareMask) << towerShareShift);
+  return packed_towerIDandShare;
+}
+
+void HGCalTriggerTowerGeometryHelper::unpackTowerIDandShare(unsigned towerIDandShare,
+                                                            int& iEta_raw,
+                                                            int& iPhi_raw,
+                                                            int& towerShare) const {
+  //eta
+  iEta_raw = (towerIDandShare >> l1t::HGCalTowerID::coord1Shift) & l1t::HGCalTowerID::coordMask;
+  unsigned iEtaSign = (towerIDandShare >> sign1Shift) & signMask;
+  iEta_raw = (iEtaSign) ? -1 * iEta_raw : iEta_raw;
+  //phi
+  iPhi_raw = (towerIDandShare >> l1t::HGCalTowerID::coord2Shift) & l1t::HGCalTowerID::coordMask;
+  unsigned iPhiSign = (towerIDandShare >> sign2Shift) & signMask;
+  iPhi_raw = (iPhiSign) ? -1 * iPhi_raw : iPhi_raw;
+  //tower share
+  towerShare = (towerIDandShare >> towerShareShift) & towerShareMask;
+}
+
+int HGCalTriggerTowerGeometryHelper::moveToCorrectSector(int iPhi_raw, int sector) const {
+  int iPhi = (iPhi_raw + sector * rotate120Deg_ + rotate180Deg_) % int(nBinsPhi_);
+  return iPhi;
+}
+
+void HGCalTriggerTowerGeometryHelper::reverseXaxis(int& iPhi) const {
+  iPhi = reverseX_ - iPhi;                          //correct x -> -x in z>0
+  iPhi = (int(nBinsPhi_) + iPhi) % int(nBinsPhi_);  // make all phi between 0 to nBinsPhi_-1
 }
 
 const std::vector<l1t::HGCalTowerCoord>& HGCalTriggerTowerGeometryHelper::getTowerCoordinates() const {
@@ -120,18 +244,76 @@ unsigned short HGCalTriggerTowerGeometryHelper::getTriggerTowerFromEtaPhi(const
   return l1t::HGCalTowerID(doNose_, zside, bin_eta, bin_phi).rawId();
 }
 
-unsigned short HGCalTriggerTowerGeometryHelper::getTriggerTower(const l1t::HGCalTriggerCell& thecell) const {
+std::unordered_map<unsigned short, float> HGCalTriggerTowerGeometryHelper::getTriggerTower(
+    const l1t::HGCalTriggerCell& thecell) const {
+  std::unordered_map<unsigned short, float> towerIDandShares = {};
   unsigned int trigger_cell_id = thecell.detId();
   // NOTE: if the TC is not found in the map than it is mapped via eta-phi coords.
   // this can be considered dangerous (silent failure of the map) but it actually allows to save
   // memory mapping explicitly only what is actually needed
   auto tower_id_itr = cells_to_trigger_towers_.find(trigger_cell_id);
-  if (tower_id_itr != cells_to_trigger_towers_.end())
-    return tower_id_itr->second;
-
-  return getTriggerTowerFromEtaPhi(thecell.position().eta(), thecell.position().phi());
+  if (tower_id_itr != cells_to_trigger_towers_.end()) {
+    towerIDandShares.insert({tower_id_itr->second, 1.0});
+    return towerIDandShares;
+  }
+  towerIDandShares.insert({getTriggerTowerFromEtaPhi(thecell.position().eta(), thecell.position().phi()), 1.0});
+  return towerIDandShares;
 }
 
-unsigned short HGCalTriggerTowerGeometryHelper::getTriggerTower(const l1t::HGCalTriggerSums& thesum) const {
-  return getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi());
+std::unordered_map<unsigned short, float> HGCalTriggerTowerGeometryHelper::getTriggerTower(
+    const l1t::HGCalTriggerSums& thesum) const {
+  std::unordered_map<unsigned short, float> towerIDandShares = {};
+  if (!splitModuleSum_) {
+    towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
+    return towerIDandShares;
+  } else {
+    HGCalTriggerModuleDetId detid(thesum.detId());
+    int moduleU = detid.moduleU();
+    int moduleV = detid.moduleV();
+    int layer = detid.layer();
+    int sector = detid.sector();
+    int zside = detid.zside();
+    int subdet = 0;
+    int splitDivisor = splitDivisorSilic_;
+    if (detid.isHScintillator()) {
+      subdet = 2;
+      splitDivisor = splitDivisorScint_;
+    } else if (detid.isEE()) {
+      subdet = 0;
+      splitDivisor = splitDivisorSilic_;
+    } else if (detid.isHSilicon()) {
+      subdet = 1;
+      splitDivisor = splitDivisorSilic_;
+    } else {  //HFNose
+      towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
+      return towerIDandShares;
+    }
+
+    unsigned packed_modID = packLayerSubdetWaferId(subdet, layer, moduleU, moduleV);
+    auto module_id_itr = modules_to_trigger_towers_.find(packed_modID);
+    if (module_id_itr != modules_to_trigger_towers_.end()) {
+      //eta variables
+      int iEta = -999;
+      int iEta_raw = -999;
+      int offsetEta = 2;
+      //phi variables
+      int iPhi = -999;
+      int iPhi_raw = -999;
+      int towerShare = -999;  //the share each tower gets from module sum
+      for (auto towerIDandShare : module_id_itr->second) {
+        unpackTowerIDandShare(towerIDandShare, iEta_raw, iPhi_raw, towerShare);
+        iEta = offsetEta + iEta_raw;
+        iPhi = moveToCorrectSector(iPhi_raw, sector);
+        if (zside == 1) {
+          reverseXaxis(iPhi);
+        }
+        towerIDandShares.insert(
+            {l1t::HGCalTowerID(doNose_, zside, iEta, iPhi).rawId(), double(towerShare) / splitDivisor});
+      }
+      return towerIDandShares;
+    } else {  // for modules not found in the mapping file (currently a few partial modules) use the traditional method.
+      towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
+      return towerIDandShares;
+    }
+  }
 }