Calo/src/HGCGuardRing.cc

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 //#define EDM_ML_DEBUG
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 }