File indexing completed on 2022-07-20 02:18:42
0001 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0004 #include "Geometry/HGCalCommonData/interface/HGCalProperty.h"
0005 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0006
0007
0008
0009 HGCalWaferType::HGCalWaferType(const std::vector<double>& rad100,
0010 const std::vector<double>& rad200,
0011 double waferSize,
0012 double zMin,
0013 int choice,
0014 unsigned int cornerCut,
0015 double cutArea)
0016 : rad100_(rad100),
0017 rad200_(rad200),
0018 waferSize_(waferSize),
0019 zMin_(zMin),
0020 choice_(choice),
0021 cutValue_(cornerCut),
0022 cutFracArea_(cutArea) {
0023 r_ = 0.5 * waferSize_;
0024 R_ = sqrt3_ * waferSize_;
0025 #ifdef EDM_ML_DEBUG
0026 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: initialized with waferR's " << waferSize_ << ":" << r_ << ":" << R_
0027 << " Choice " << choice_ << " Cuts " << cutValue_ << ":" << cutFracArea_ << " zMin "
0028 << zMin_ << " with " << rad100_.size() << ":" << rad200_.size() << " parameters for R:";
0029 for (unsigned k = 0; k < rad100_.size(); ++k)
0030 edm::LogVerbatim("HGCalGeom") << "[" << k << "] 100:200 " << rad100_[k] << " 200:300 " << rad200_[k];
0031 #endif
0032 }
0033
0034 int HGCalWaferType::getCassette(int index, const HGCalParameters::waferInfo_map& wafers) {
0035 auto itr = wafers.find(index);
0036 return ((itr == wafers.end()) ? -1 : ((itr->second).cassette));
0037 }
0038
0039 int HGCalWaferType::getOrient(int index, const HGCalParameters::waferInfo_map& wafers) {
0040 auto itr = wafers.find(index);
0041 return ((itr == wafers.end()) ? -1 : ((itr->second).orient));
0042 }
0043
0044 int HGCalWaferType::getPartial(int index, const HGCalParameters::waferInfo_map& wafers) {
0045 auto itr = wafers.find(index);
0046 return ((itr == wafers.end()) ? -1 : ((itr->second).part));
0047 }
0048
0049 int HGCalWaferType::getType(int index, const HGCalParameters::waferInfo_map& wafers) {
0050 auto itr = wafers.find(index);
0051 return ((itr == wafers.end()) ? -1 : ((itr->second).type));
0052 }
0053
0054 int HGCalWaferType::getType(int index, const std::vector<int>& indices, const std::vector<int>& properties) {
0055 auto itr = std::find(std::begin(indices), std::end(indices), index);
0056 int type = (itr == std::end(indices))
0057 ? -1
0058 : HGCalProperty::waferThick(properties[static_cast<unsigned int>(itr - std::begin(indices))]);
0059 return type;
0060 }
0061
0062 int HGCalWaferType::getType(double xpos, double ypos, double zpos) {
0063 std::vector<double> xc(HGCalParameters::k_CornerSize, 0);
0064 std::vector<double> yc(HGCalParameters::k_CornerSize, 0);
0065 xc[0] = xpos + r_;
0066 yc[0] = ypos + 0.5 * R_;
0067 xc[1] = xpos;
0068 yc[1] = ypos + R_;
0069 xc[2] = xpos - r_;
0070 yc[2] = ypos + 0.5 * R_;
0071 xc[3] = xpos - r_;
0072 yc[3] = ypos - 0.5 * R_;
0073 xc[4] = xpos;
0074 yc[4] = ypos - R_;
0075 xc[5] = xpos + r_;
0076 yc[5] = ypos - 0.5 * R_;
0077 const auto& rv = rLimits(zpos);
0078 std::vector<int> fine, coarse;
0079 for (unsigned int k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0080 double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0081 if (rpos <= rv.first)
0082 fine.emplace_back(k);
0083 else if (rpos <= rv.second)
0084 coarse.emplace_back(k);
0085 }
0086 int type(-2);
0087 double fracArea(0);
0088 if (choice_ == 1) {
0089 if (fine.size() >= cutValue_)
0090 type = HGCSiliconDetId::HGCalFine;
0091 else if (coarse.size() >= cutValue_)
0092 type = HGCSiliconDetId::HGCalCoarseThin;
0093 else
0094 type = HGCSiliconDetId::HGCalCoarseThick;
0095 } else {
0096 if (fine.size() >= 4)
0097 type = HGCSiliconDetId::HGCalFine;
0098 else if (coarse.size() >= 4 && fine.size() <= 1)
0099 type = HGCSiliconDetId::HGCalCoarseThin;
0100 else if (coarse.size() < 2 && fine.empty())
0101 type = HGCSiliconDetId::HGCalCoarseThick;
0102 else if (!fine.empty())
0103 type = -1;
0104 if (type <= -1) {
0105 unsigned int kmax = (type == -1) ? fine.size() : coarse.size();
0106 std::vector<double> xcn, ycn;
0107 for (unsigned int k = 0; k < kmax; ++k) {
0108 unsigned int k1 = (type == -1) ? fine[k] : coarse[k];
0109 unsigned int k2 = (k1 == xc.size() - 1) ? 0 : k1 + 1;
0110 bool ok = ((type == -1) ? (std::find(fine.begin(), fine.end(), k2) != fine.end())
0111 : (std::find(coarse.begin(), coarse.end(), k2) != coarse.end()));
0112 xcn.emplace_back(xc[k1]);
0113 ycn.emplace_back(yc[k1]);
0114 if (!ok) {
0115 double rr = (type == -1) ? rv.first : rv.second;
0116 const auto& xy = intersection(k1, k2, xc, yc, xpos, ypos, rr);
0117 xcn.emplace_back(xy.first);
0118 ycn.emplace_back(xy.second);
0119 }
0120 }
0121 fracArea = areaPolygon(xcn, ycn) / areaPolygon(xc, yc);
0122 type = (fracArea > cutFracArea_) ? -(1 + type) : -type;
0123 }
0124 }
0125 #ifdef EDM_ML_DEBUG
0126 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: position " << xpos << ":" << ypos << ":" << zpos << " R "
0127 << ":" << rv.first << ":" << rv.second << " corners|type " << fine.size() << ":"
0128 << coarse.size() << ":" << fracArea << ":" << type;
0129 #endif
0130 return type;
0131 }
0132
0133 std::pair<double, double> HGCalWaferType::rLimits(double zpos) {
0134 double zz = std::abs(zpos);
0135 if (zz < zMin_)
0136 zz = zMin_;
0137 zz *= HGCalParameters::k_ScaleFromDDD;
0138 double rfine = rad100_[0];
0139 double rcoarse = rad200_[0];
0140 for (int i = 1; i < 5; ++i) {
0141 rfine *= zz;
0142 rfine += rad100_[i];
0143 rcoarse *= zz;
0144 rcoarse += rad200_[i];
0145 }
0146 rfine *= HGCalParameters::k_ScaleToDDD;
0147 rcoarse *= HGCalParameters::k_ScaleToDDD;
0148 #ifdef EDM_ML_DEBUG
0149 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: Z " << zpos << ":" << zz << " R " << rfine << ":" << rcoarse;
0150 #endif
0151 return std::make_pair(rfine, rcoarse);
0152 }
0153
0154 double HGCalWaferType::areaPolygon(std::vector<double> const& x, std::vector<double> const& y) {
0155 double area = 0.0;
0156 int n = static_cast<int>(x.size());
0157 int j = n - 1;
0158 for (int i = 0; i < n; ++i) {
0159 area += ((x[j] + x[i]) * (y[i] - y[j]));
0160 j = i;
0161 }
0162 return (0.5 * area);
0163 }
0164
0165 std::pair<double, double> HGCalWaferType::intersection(
0166 int k1, int k2, std::vector<double> const& x, std::vector<double> const& y, double xpos, double ypos, double rr) {
0167 double slope = (x[k1] - x[k2]) / (y[k1] - y[k2]);
0168 double interc = x[k1] - slope * y[k1];
0169 double xx[2], yy[2], dist[2];
0170 double v1 = std::sqrt((slope * slope + 1) * rr * rr - (interc * interc));
0171 yy[0] = (-slope * interc + v1) / (1 + slope * slope);
0172 yy[1] = (-slope * interc - v1) / (1 + slope * slope);
0173 for (int i = 0; i < 2; ++i) {
0174 xx[i] = (slope * yy[i] + interc);
0175 dist[i] = ((xx[i] - xpos) * (xx[i] - xpos)) + ((yy[i] - ypos) * (yy[i] - ypos));
0176 }
0177 #ifdef EDM_ML_DEBUG
0178 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: InterSection " << dist[0] << ":" << xx[0] << ":" << yy[0] << " vs "
0179 << dist[1] << ":" << xx[1] << ":" << yy[1];
0180 #endif
0181 if (dist[0] > dist[1])
0182 return std::make_pair(xx[1], yy[1]);
0183 else
0184 return std::make_pair(xx[0], yy[0]);
0185 }
