Merge pull request #41380 from bsunanda/Phase2-hgx337B

Phase2-hgx337B Try to solve the cassette shift issue and fix a bug in Geometry/HGCalCommonData

Geometry/HGCalCommonData/interface/HGCalCassette.h

@@ -13,10 +13,11 @@ class HGCalCassette {
   void setParameter(int cassette, const std::vector<double>& shifts);
   std::pair<double, double> getShift(int layer, int zside, int cassette) const;
   static int cassetteIndex(int det, int layer, int zside, int cassette);
+  static int cassetteType(int det, int zside, int cassette);
 
 private:
-  const std::vector<int> positEE_ = {2, 1, 0, 5, 4, 3};
-  const std::vector<int> positHE_ = {5, 4, 3, 2, 1, 0, 11, 10, 9, 8, 7, 6};
+  static constexpr int positEE_[6] = {2, 1, 0, 5, 4, 3};
+  static constexpr int positHE_[12] = {5, 4, 3, 2, 1, 0, 11, 10, 9, 8, 7, 6};
   int cassette_;
   bool typeHE_;
   std::vector<double> shifts_;

Geometry/HGCalCommonData/plugins/DDHGCalMixRotatedLayer.cc

@@ -469,12 +469,13 @@ void DDHGCalMixRotatedLayer::positionMix(const DDLogicalPart& glog,
       double phi2 = dphi * (fimax - fimin + 1);
       auto cshift = cassette_.getShift(layer + 1, 1, cassette);
 #ifdef EDM_ML_DEBUG
+      int cassette0 = HGCalCassette::cassetteType(2, 1, cassette);  //
       edm::LogVerbatim("HGCalGeom") << "DDHGCalMixRotatedLayer: Layer " << copy << ":" << (layer + 1) << " iR "
                                     << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
                                     << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R " << r1 << ":"
                                     << r2 << " Thick " << (2.0 * hthickl) << " phi " << fimin << ":" << fimax << ":"
                                     << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2) << " cassette " << cassette
-                                    << " Shift " << cshift.first << ":" << cshift.second;
+                                    << ":" << cassette0 << " Shift " << cshift.first << ":" << cshift.second;
 #endif
       std::string name = namesTop_[ii] + "L" + std::to_string(copy) + "F" + std::to_string(k);
       ++k;

Geometry/HGCalCommonData/plugins/DDHGCalSiliconRotatedModule.cc

@@ -427,10 +427,10 @@ void DDHGCalSiliconRotatedModule::positionSensitive(const DDLogicalPart& glog, i
     double xorig = xyoff.first + nc * delx;
     double yorig = xyoff.second + nr * dy;
     double angle = std::atan2(yorig, xorig);
-    edm::LogVerbatim("HGCalGeom") << "DDHGCalSiliconRotatedLayer::Wafer: layer " << layer + 1 << " cassette "
+    edm::LogVerbatim("HGCalGeom") << "DDHGCalSiliconRotatedModule::Wafer: layer " << layer + 1 << " cassette "
                                   << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original "
                                   << xorig << ":" << yorig << ":" << convertRadToDeg(angle) << " Final " << xpos << ":"
-                                  << ypos << " u|v " << u << ":" << v << " typee|part|orient|place " << type << ":"
+                                  << ypos << " u|v " << u << ":" << v << " type|part|orient|place " << type << ":"
                                   << part << ":" << orien << ":" << place;
 #endif
     std::string wafer;

Geometry/HGCalCommonData/plugins/dd4hep/DDHGCalMixRotatedLayer.cc

@@ -396,13 +396,15 @@ struct HGCalMixRotatedLayer {
         double phi2 = (forFireworks_ == 1) ? (dphi * (fimax - fimin + 1)) : (dphi * fimax);
         auto cshift = cassette_.getShift(layer + 1, 1, cassette);
 #ifdef EDM_ML_DEBUG
+        int cassette0 = HGCalCassette::cassetteType(2, 1, cassette);  //
         edm::LogVerbatim("HGCalGeom") << "DDHGCalMixRotatedLayer: Layer " << copy << ":" << (layer + 1) << " iR "
                                       << std::get<1>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << ":"
                                       << std::get<2>(HGCalTileIndex::tileUnpack(tileIndex_[ti])) << " R "
                                       << cms::convert2mm(r1) << ":" << cms::convert2mm(r2) << " Thick "
                                       << cms::convert2mm((2.0 * hthickl)) << " phi " << fimin << ":" << fimax << ":"
                                       << convertRadToDeg(phi1) << ":" << convertRadToDeg(phi2) << " cassette "
-                                      << cassette << " Shift " << cshift.first << ":" << cshift.second;
+                                      << cassette << ":" << cassett0 << " Shift " << cshift.first << ":"
+                                      << cshift.second;
 #endif
         std::string name = namesTop_[ii] + "L" + std::to_string(copy) + "F" + std::to_string(k);
         ++k;

Geometry/HGCalCommonData/plugins/dd4hep/DDHGCalSiliconRotatedModule.cc

@@ -367,10 +367,10 @@ struct HGCalSiliconRotatedModule {
       double xorig = xyoff.first + nc * delx;
       double yorig = xyoff.second + nr * dy;
       double angle = std::atan2(yorig, xorig);
-      edm::LogVerbatim("HGCalGeom") << "DDHGCalSiliconRotatedLayer::Wafer: layer " << layer + 1 << " cassette "
+      edm::LogVerbatim("HGCalGeom") << "DDHGCalSiliconRotatedModule::Wafer: layer " << layer + 1 << " cassette "
                                     << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original "
                                     << xorig << ":" << yorig << ":" << convertRadToDeg(angle) << " Final " << xpos
-                                    << ":" << ypos << " u|v " << u << ":" << v << " typee|part|orient|place " << type
+                                    << ":" << ypos << " u|v " << u << ":" << v << " type|part|orient|place " << type
                                     << ":" << part << ":" << orien << ":" << place;
 #endif
       std::string wafer;

Geometry/HGCalCommonData/src/HGCalCassette.cc

@@ -28,11 +28,11 @@ void HGCalCassette::setParameter(int cassette, const std::vector<double>& shifts
 std::pair<double, double> HGCalCassette::getShift(int layer, int zside, int cassette) const {
   int locc = (zside < 0) ? (cassette - 1) : (typeHE_ ? positHE_[cassette - 1] : positEE_[cassette - 1]);
   int loc = 2 * (cassette_ * (layer - 1) + locc);
-  std::pair<double, double> xy = std::make_pair(shifts_[loc], shifts_[loc + 1]);
+  std::pair<double, double> xy = std::make_pair(-zside * shifts_[loc], shifts_[loc + 1]);
 #ifdef EDM_ML_DEBUG
-  edm::LogVerbatim("HGCalGeom") << "HGCalCassette::getShift: Layer " << layer << " zside " << zside << " cassette "
-                                << cassette << " Loc " << locc << ":" << loc << " shift " << xy.first << ":"
-                                << xy.second;
+  edm::LogVerbatim("HGCalGeom") << "HGCalCassette::getShift: Layer " << layer << " zside " << zside << " type "
+                                << typeHE_ << " cassette " << cassette << " Loc " << locc << ":" << loc << " shift "
+                                << xy.first << ":" << xy.second;
 #endif
   return xy;
 }
@@ -41,3 +41,8 @@ int HGCalCassette::cassetteIndex(int det, int layer, int side, int cassette) {
   int zs = (side > 0) ? factor_ : 0;
   return (((zs + det) * factor_ + layer) * factor_ + cassette);
 }
+
+int HGCalCassette::cassetteType(int det, int zside, int cassette) {
+  int type = (zside < 0) ? cassette : ((det == 0) ? (1 + positEE_[cassette - 1]) : (1 + positHE_[cassette - 1]));
+  return type;
+}

Geometry/HGCalCommonData/src/HGCalDDDConstants.cc

@@ -189,24 +189,28 @@ std::array<int, 3> HGCalDDDConstants::assignCellTrap(float x, float y, float z,
   int ll = layer - hgpar_->firstLayer_;
   xx -= hgpar_->xLayerHex_[ll];
   yy -= hgpar_->yLayerHex_[ll];
-  if (mode_ == HGCalGeometryMode::TrapezoidCassette) {
-    int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
-    auto cshift = hgcassette_.getShift(layer, zside, cassette);
-    int invert = ((layerType(layer) == HGCalTypes::WaferCenter) && (zside > 0)) ? 1 : -1;
-#ifdef EDM_ML_DEBUG
-    edm::LogVerbatim("HGCalGeomT") << "Cassette " << cassette << " Shift " << (invert * -cshift.first) << ":"
-                                   << cshift.second;
-#endif
-    xx -= (invert * cshift.first);
-    yy -= cshift.second;
-  }
   double phi = (((yy == 0.0) && (xx == 0.0)) ? 0. : std::atan2(yy, xx));
   if (phi < 0)
     phi += (2.0 * M_PI);
   if (indx.second != 0)
     iphi = (1 + static_cast<int>(phi / indx.second)) % hgpar_->scintCells(layer);
   if (iphi == 0)
     iphi = hgpar_->scintCells(layer);
+  if (mode_ == HGCalGeometryMode::TrapezoidCassette) {
+    int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
+    auto cshift = hgcassette_.getShift(layer, -1, cassette);
+#ifdef EDM_ML_DEBUG
+    std::ostringstream st1;
+    st1 << "Cassette " << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original " << xx << ":"
+        << yy;
+#endif
+    xx += (zside * cshift.first);
+    yy -= cshift.second;
+#ifdef EDM_ML_DEBUG
+    st1 << " Shifted " << xx << ":" << yy;
+    edm::LogVerbatim("HGCalGeomT") << st1.str();
+#endif
+  }
   type = hgpar_->scintType(layer);
   double r = std::sqrt(xx * xx + yy * yy);
   auto ir = std::lower_bound(hgpar_->radiusLayer_[type].begin(), hgpar_->radiusLayer_[type].end(), r);
@@ -753,18 +757,22 @@ std::pair<float, float> HGCalDDDConstants::localToGlobal8(
   int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
   auto ktr = hgpar_->waferInfoMap_.find(indx);
   if ((mode_ == HGCalGeometryMode::Hexagon8Cassette) && (ktr != hgpar_->waferInfoMap_.end())) {
-    int invert = ((layerType(lay) == HGCalTypes::WaferCenter) && (zside > 0)) ? 1 : -1;
-    auto cshift = hgcassette_.getShift(lay, zside, (ktr->second).cassette);
+    auto cshift = hgcassette_.getShift(lay, -1, (ktr->second).cassette);
+    std::ostringstream st1;
     if (debug)
-      edm::LogVerbatim("HGCalGeom") << "Cassette " << (ktr->second).cassette << " Shift " << (invert * cshift.first)
-                                    << ":" << cshift.second;
+      st1 << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second << " Original "
+          << x << ":" << y;
     if (!reco) {
-      x += ((HGCalParameters::k_ScaleToDDD)*invert * cshift.first);
+      x -= ((HGCalParameters::k_ScaleToDDD)*zside * cshift.first);
       y += ((HGCalParameters::k_ScaleToDDD)*cshift.second);
     } else {
-      x += (invert * cshift.first);
+      x -= (zside * cshift.first);
       y += cshift.second;
     }
+    if (debug) {
+      st1 << " Final " << x << ":" << y;
+      edm::LogVerbatim("HGCalGeom") << st1.str();
+    }
   }
   if (debug)
     edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
@@ -862,18 +870,22 @@ std::pair<float, float> HGCalDDDConstants::locateCell(
     x += xy.first;
     y += xy.second;
     if ((mode_ == HGCalGeometryMode::Hexagon8Cassette) && (ktr != hgpar_->waferInfoMap_.end())) {
-      int invert = ((layerType(lay) == HGCalTypes::WaferCenter) && (zside > 0)) ? 1 : -1;
-      auto cshift = hgcassette_.getShift(lay, zside, (ktr->second).cassette);
+      auto cshift = hgcassette_.getShift(lay, -1, (ktr->second).cassette);
+      std::ostringstream st1;
       if (debug)
-        edm::LogVerbatim("HGCalGeom") << "Cassette " << (ktr->second).cassette << " Shift " << (invert * cshift.first)
-                                      << ":" << cshift.second;
+        st1 << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second
+            << " Original " << x << ":" << y;
       if (!reco) {
-        x += ((HGCalParameters::k_ScaleToDDD)*invert * cshift.first);
+        x -= ((HGCalParameters::k_ScaleToDDD)*cshift.first);
         y += ((HGCalParameters::k_ScaleToDDD)*cshift.second);
       } else {
-        x += (invert * cshift.first);
+        x -= cshift.first;
         y += cshift.second;
       }
+      if (debug) {
+        st1 << " Final " << x << ":" << y;
+        edm::LogVerbatim("HGCalGeom") << st1.str();
+      }
     }
     if (debug)
       edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
@@ -930,18 +942,22 @@ std::pair<float, float> HGCalDDDConstants::locateCellTrap(
     int ll = lay - hgpar_->firstLayer_;
     x += hgpar_->xLayerHex_[ll];
     y += hgpar_->yLayerHex_[ll];
+    if (irad < 0)
+      x = -x;
     if (mode_ == HGCalGeometryMode::TrapezoidCassette) {
       int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
-      int invert = ((layerType(lay) == HGCalTypes::WaferCenter) && (zside > 0)) ? 1 : -1;
-      auto cshift = hgcassette_.getShift(lay, zside, cassette);
+      auto cshift = hgcassette_.getShift(lay, -1, cassette);
+      std::ostringstream st1;
       if (debug)
-        edm::LogVerbatim("HGCalGeom") << "Cassette " << cassette << " Shift " << (invert * cshift.first) << ":"
-                                      << cshift.second;
-      x += (invert * cshift.first);
+        st1 << "Cassette " << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original " << x << ":"
+            << y;
+      x -= cshift.first;
       y += cshift.second;
+      if (debug) {
+        st1 << " Final " << x << ":" << y;
+        edm::LogVerbatim("HGCalGeom") << st1.str();
+      }
     }
-    if (irad < 0)
-      x = -x;
   }
   if (!reco) {
     x *= HGCalParameters::k_ScaleToDDD;
@@ -1465,12 +1481,11 @@ void HGCalDDDConstants::waferFromPosition(const double x,
       int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
       auto ktr = hgpar_->waferInfoMap_.find(indx);
       if (ktr != hgpar_->waferInfoMap_.end()) {
-        int invert = ((layerType(layer) == HGCalTypes::WaferCenter) && (zside > 0)) ? 1 : -1;
-        auto cshift = hgcassette_.getShift(layer, zside, (ktr->second).cassette);
+        auto cshift = hgcassette_.getShift(layer, -1, (ktr->second).cassette);
         if (debug)
-          edm::LogVerbatim("HGCalGeom") << "Cassette " << (ktr->second).cassette << " Shift " << (invert * cshift.first)
-                                        << ":" << cshift.second;
-        dx0 = (invert * cshift.first);
+          edm::LogVerbatim("HGCalGeom") << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":"
+                                        << cshift.second;
+        dx0 = -cshift.first;
         dy0 = cshift.second;
       }
     }

Geometry/HGCalCommonData/test/HGCalParameterTester.cc

@@ -134,8 +134,8 @@ void HGCalParameterTester::analyze(const edm::Event& iEvent, const edm::EventSet
     edm::LogVerbatim("HGCalGeom") << "UseSimWt: " << phgp->useSimWt_;
     edm::LogVerbatim("HGCalGeom") << "Wafer Parameters: " << phgp->waferSize_ << ":" << phgp->waferR_ << ":"
                                   << phgp->waferThick_ << ":" << phgp->sensorSeparation_ << ":"
-                                  << phgp->sensorSizeOffset_ << ":" << phgp->guardRingOffset_ << ":"
-                                  << phgp->mouseBite_;
+                                  << phgp->sensorSizeOffset_ << ":" << phgp->guardRingOffset_ << ":" << phgp->mouseBite_
+                                  << ":" << phgp->useOffset_;
     myPrint("waferThickness", phgp->waferThickness_, 10);
     edm::LogVerbatim("HGCalGeom") << "nCells_: " << phgp->nCellsFine_ << ":" << phgp->nCellsCoarse_;
     edm::LogVerbatim("HGCalGeom") << "nSectors_: " << phgp->nSectors_;
@@ -145,6 +145,7 @@ void HGCalParameterTester::analyze(const edm::Event& iEvent, const edm::EventSet
     edm::LogVerbatim("HGCalGeom") << "mode_: " << phgp->mode_;
     edm::LogVerbatim("HGCalGeom") << "cassettes_: " << phgp->cassettes_;
 
+    edm::LogVerbatim("HGCalGeom") << "waferUVMax: " << phgp->waferUVMax_;
     myPrint("waferUVMaxLayer", phgp->waferUVMaxLayer_, 20);
     myPrint("CellThickness", phgp->cellThickness_, 10);
     myPrint("radius100to200", phgp->radius100to200_, 10);

Geometry/HGCalCommonData/test/python/testHGCalParametersDDD_cfg.py

@@ -2,7 +2,7 @@
 # Way to use this:
 #   cmsRun testHGCalParametersDDD_cfg.py type=V17
 #
-#   Options for type V16, V17, V17n
+#   Options for type V16, V17, V17n, V17Shift
 #
 ###############################################################################
 import FWCore.ParameterSet.Config as cms
@@ -16,7 +16,7 @@
                  "V17",
                   VarParsing.VarParsing.multiplicity.singleton,
                   VarParsing.VarParsing.varType.string,
-                  "type of operations: V16, V17, V17n")
+                  "type of operations: V16, V17, V17n, V17Shift")
 
 ### get and parse the command line arguments
 options.parseArguments()

SimG4CMS/Calo/src/HGCalSD.cc

@@ -210,25 +210,27 @@ uint32_t HGCalSD::setDetUnitId(const G4Step* aStep) {
     double xx = (hid1.zside() > 0) ? xy.first : -xy.first;
     double dx = xx - (hitPoint.x() / CLHEP::cm);
     double dy = xy.second - (hitPoint.y() / CLHEP::cm);
-    double diff = std::sqrt(dx * dx + dy * dy);
-    constexpr double tol = 2.0;
+    double diff = (dx * dx + dy * dy);
+    constexpr double tol = 2.0 * 2.0;
     bool valid1 = hgcons_->isValidHex8(hid1.layer(), hid1.waferU(), hid1.waferV(), hid1.cellU(), hid1.cellV(), true);
     if ((diff > tol) || (!valid1))
       pid = "HGCalError";
     auto partn = hgcons_->waferTypeRotation(hid1.layer(), hid1.waferU(), hid1.waferV(), false, false);
     int indx = HGCalWaferIndex::waferIndex(layer, hid1.waferU(), hid1.waferV());
-    double phi = std::atan2(hitPoint.y(), hitPoint.x());
     edm::LogVerbatim(pid) << "CheckID " << HGCSiliconDetId(id) << " Layer:Module:Cell:ModuleLev " << layer << ":"
                           << module << ":" << cell << ":" << moduleLev << " SimWt:history " << useSimWt_ << ":"
                           << touch->GetHistoryDepth() << ":" << levelT1_ << ":" << levelT2_ << " input position: ("
-                          << hitPoint.x() / CLHEP::cm << ", " << hitPoint.y() / CLHEP::cm << ":" << convertRadToDeg(phi)
-                          << "); position from ID (" << xx << ", " << xy.second << ") distance " << dx << ":" << dy
-                          << ":" << diff << " Valid " << valid1 << " Wafer type|rotation " << partn.first << ":"
-                          << partn.second << " Part:Orient:Cassette " << std::get<1>(hgcons_->waferFileInfo(indx))
-                          << ":" << std::get<2>(hgcons_->waferFileInfo(indx)) << ":"
+                          << hitPoint.x() / CLHEP::cm << ", " << hitPoint.y() / CLHEP::cm << ":"
+                          << convertRadToDeg(std::atan2(hitPoint.y(), hitPoint.x())) << "); position from ID (" << xx
+                          << ", " << xy.second << ") distance " << dx << ":" << dy << ":" << std::sqrt(diff)
+                          << " Valid " << valid1 << " Wafer type|rotation " << partn.first << ":" << partn.second
+                          << " Part:Orient:Cassette " << std::get<1>(hgcons_->waferFileInfo(indx)) << ":"
+                          << std::get<2>(hgcons_->waferFileInfo(indx)) << ":"
                           << std::get<3>(hgcons_->waferFileInfo(indx)) << " CassetteShift " << cshift;
-    xy = hgcons_->locateCell(hid1, true);
-    printDetectorLevels(touch);
+    if ((diff > tol) || (!valid1)) {
+      printDetectorLevels(touch);
+      hgcons_->locateCell(hid1, true);
+    }
   }
   return id;
 }