File indexing completed on 2023-05-25 00:54:20
0001 #include "SimG4CMS/Calo/interface/HGCGuardRing.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0004 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0005 #include <iostream>
0006
0007
0008
0009 HGCGuardRing::HGCGuardRing(const HGCalDDDConstants& hgc)
0010 : hgcons_(hgc),
0011 modeUV_(hgcons_.geomMode()),
0012 waferSize_(hgcons_.waferSize(false)),
0013 sensorSizeOffset_(hgcons_.getParameter()->sensorSizeOffset_),
0014 guardRingOffset_(hgcons_.getParameter()->guardRingOffset_) {
0015 offset_ = sensorSizeOffset_ + 2.0 * guardRingOffset_;
0016 xmax_ = 0.5 * (waferSize_ - offset_);
0017 ymax_ = xmax_ / sqrt3_;
0018 #ifdef EDM_ML_DEBUG
0019 edm::LogVerbatim("HGCSim") << "Creating HGCGuardRing with wafer size " << waferSize_ << ", Offsets "
0020 << sensorSizeOffset_ << ":" << guardRingOffset_ << ":" << offset_ << ", and mode "
0021 << modeUV_ << " xmax|ymax " << xmax_ << ":" << ymax_;
0022 #endif
0023 }
0024
0025 bool HGCGuardRing::exclude(G4ThreeVector& point, int zside, int frontBack, int layer, int waferU, int waferV) {
0026 bool check(false);
0027 if ((modeUV_ == HGCalGeometryMode::Hexagon8Module) || (modeUV_ == HGCalGeometryMode::Hexagon8Cassette)) {
0028 int index = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
0029 int partial = HGCalWaferType::getPartial(index, hgcons_.getParameter()->waferInfoMap_);
0030 if (partial == HGCalTypes::WaferFull) {
0031 double dx = std::abs(point.x());
0032 double dy = std::abs(point.y());
0033 if (dx > xmax_) {
0034 check = true;
0035 } else if (dy > (2 * ymax_)) {
0036 check = true;
0037 } else {
0038 check = (dx > (sqrt3_ * (2 * ymax_ - dy)));
0039 }
0040 #ifdef EDM_ML_DEBUG
0041 edm::LogVerbatim("HGCSim") << "HGCGuardRing:: Point " << point << " zside " << zside << " layer " << layer
0042 << " wafer " << waferU << ":" << waferV << " partial type " << partial << ":"
0043 << HGCalTypes::WaferFull << " x " << dx << ":" << xmax_ << " y " << dy << ":" << ymax_
0044 << " check " << check;
0045 #endif
0046 } else {
0047 int orient = HGCalWaferType::getOrient(index, hgcons_.getParameter()->waferInfoMap_);
0048 if (modeUV_ == HGCalGeometryMode::Hexagon8Module) {
0049 std::vector<std::pair<double, double> > wxy =
0050 HGCalWaferMask::waferXY(partial, orient, zside, waferSize_, offset_, 0.0, 0.0);
0051 check = !(insidePolygon(point.x(), point.y(), wxy));
0052 #ifdef EDM_ML_DEBUG
0053 std::ostringstream st1;
0054 st1 << "HGCGuardRing:: Point " << point << " Partial/orient/zside/size/offset " << partial << ":" << orient
0055 << ":" << zside << ":" << waferSize_ << offset_ << " with " << wxy.size() << " points:";
0056 for (unsigned int k = 0; k < wxy.size(); ++k)
0057 st1 << " (" << wxy[k].first << ", " << wxy[k].second << ")";
0058 edm::LogVerbatim("HGCSim") << st1.str();
0059 #endif
0060 } else {
0061 int placement = HGCalCell::cellPlacementIndex(zside, frontBack, orient);
0062 std::vector<std::pair<double, double> > wxy =
0063 HGCalWaferMask::waferXY(partial, placement, waferSize_, offset_, 0.0, 0.0);
0064 check = !(insidePolygon(point.x(), point.y(), wxy));
0065 #ifdef EDM_ML_DEBUG
0066 std::ostringstream st1;
0067 st1 << "HGCGuardRing:: Point " << point << " Partial/frontback/orient/zside/placeemnt/size/offset " << partial
0068 << ":" << frontBack << ":" << orient << ":" << zside << ":" << placement << ":" << waferSize_ << offset_
0069 << " with " << wxy.size() << " points:";
0070 for (unsigned int k = 0; k < wxy.size(); ++k)
0071 st1 << " (" << wxy[k].first << ", " << wxy[k].second << ")";
0072 edm::LogVerbatim("HGCSim") << st1.str();
0073 #endif
0074 }
0075 }
0076 }
0077 return check;
0078 }
0079
0080 bool HGCGuardRing::insidePolygon(double x, double y, const std::vector<std::pair<double, double> >& xyv) {
0081 int counter(0);
0082 double x1(xyv[0].first), y1(xyv[0].second);
0083 for (unsigned i1 = 1; i1 <= xyv.size(); i1++) {
0084 unsigned i2 = (i1 % xyv.size());
0085 double x2(xyv[i2].first), y2(xyv[i2].second);
0086 if (y > std::min(y1, y2)) {
0087 if (y <= std::max(y1, y2)) {
0088 if (x <= std::max(x1, x2)) {
0089 if (y1 != y2) {
0090 double xinter = (y - y1) * (x2 - x1) / (y2 - y1) + x1;
0091 if ((x1 == x2) || (x <= xinter))
0092 ++counter;
0093 }
0094 }
0095 }
0096 }
0097 x1 = x2;
0098 y1 = y2;
0099 }
0100
0101 if (counter % 2 == 0)
0102 return false;
0103 else
0104 return true;
0105 }