File indexing completed on 2021-02-14 12:48:27
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 HGCalWaferType::~HGCalWaferType() {}
0035
0036 int HGCalWaferType::getType(double xpos, double ypos, double zpos) {
0037 std::vector<double> xc(HGCalParameters::k_CornerSize, 0);
0038 std::vector<double> yc(HGCalParameters::k_CornerSize, 0);
0039 xc[0] = xpos + r_;
0040 yc[0] = ypos + 0.5 * R_;
0041 xc[1] = xpos;
0042 yc[1] = ypos + R_;
0043 xc[2] = xpos - r_;
0044 yc[2] = ypos + 0.5 * R_;
0045 xc[3] = xpos - r_;
0046 yc[3] = ypos - 0.5 * R_;
0047 xc[4] = xpos;
0048 yc[4] = ypos - R_;
0049 xc[5] = xpos + r_;
0050 yc[5] = ypos - 0.5 * R_;
0051 const auto& rv = rLimits(zpos);
0052 std::vector<int> fine, coarse;
0053 for (unsigned int k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0054 double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0055 if (rpos <= rv.first)
0056 fine.emplace_back(k);
0057 else if (rpos <= rv.second)
0058 coarse.emplace_back(k);
0059 }
0060 int type(-2);
0061 double fracArea(0);
0062 if (choice_ == 1) {
0063 if (fine.size() >= cutValue_)
0064 type = HGCSiliconDetId::HGCalFine;
0065 else if (coarse.size() >= cutValue_)
0066 type = HGCSiliconDetId::HGCalCoarseThin;
0067 else
0068 type = HGCSiliconDetId::HGCalCoarseThick;
0069 } else {
0070 if (fine.size() >= 4)
0071 type = HGCSiliconDetId::HGCalFine;
0072 else if (coarse.size() >= 4 && fine.size() <= 1)
0073 type = HGCSiliconDetId::HGCalCoarseThin;
0074 else if (coarse.size() < 2 && fine.empty())
0075 type = HGCSiliconDetId::HGCalCoarseThick;
0076 else if (!fine.empty())
0077 type = -1;
0078 if (type <= -1) {
0079 unsigned int kmax = (type == -1) ? fine.size() : coarse.size();
0080 std::vector<double> xcn, ycn;
0081 for (unsigned int k = 0; k < kmax; ++k) {
0082 unsigned int k1 = (type == -1) ? fine[k] : coarse[k];
0083 unsigned int k2 = (k1 == xc.size() - 1) ? 0 : k1 + 1;
0084 bool ok = ((type == -1) ? (std::find(fine.begin(), fine.end(), k2) != fine.end())
0085 : (std::find(coarse.begin(), coarse.end(), k2) != coarse.end()));
0086 xcn.emplace_back(xc[k1]);
0087 ycn.emplace_back(yc[k1]);
0088 if (!ok) {
0089 double rr = (type == -1) ? rv.first : rv.second;
0090 const auto& xy = intersection(k1, k2, xc, yc, xpos, ypos, rr);
0091 xcn.emplace_back(xy.first);
0092 ycn.emplace_back(xy.second);
0093 }
0094 }
0095 fracArea = areaPolygon(xcn, ycn) / areaPolygon(xc, yc);
0096 type = (fracArea > cutFracArea_) ? -(1 + type) : -type;
0097 }
0098 }
0099 #ifdef EDM_ML_DEBUG
0100 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: position " << xpos << ":" << ypos << ":" << zpos << " R "
0101 << ":" << rv.first << ":" << rv.second << " corners|type " << fine.size() << ":"
0102 << coarse.size() << ":" << fracArea << ":" << type;
0103 #endif
0104 return type;
0105 }
0106
0107 int HGCalWaferType::getType(int index, const std::vector<int>& indices, const std::vector<int>& properties) {
0108 auto itr = std::find(std::begin(indices), std::end(indices), index);
0109 int type = (itr == std::end(indices))
0110 ? -1
0111 : HGCalProperty::waferThick(properties[static_cast<unsigned int>(itr - std::begin(indices))]);
0112 return type;
0113 }
0114
0115 int HGCalWaferType::getType(int index, const HGCalParameters::waferInfo_map& wafers) {
0116 auto itr = wafers.find(index);
0117 return ((itr == wafers.end()) ? -1 : ((itr->second).type));
0118 }
0119
0120 std::pair<double, double> HGCalWaferType::rLimits(double zpos) {
0121 double zz = std::abs(zpos);
0122 if (zz < zMin_)
0123 zz = zMin_;
0124 zz *= HGCalParameters::k_ScaleFromDDD;
0125 double rfine = rad100_[0];
0126 double rcoarse = rad200_[0];
0127 for (int i = 1; i < 5; ++i) {
0128 rfine *= zz;
0129 rfine += rad100_[i];
0130 rcoarse *= zz;
0131 rcoarse += rad200_[i];
0132 }
0133 rfine *= HGCalParameters::k_ScaleToDDD;
0134 rcoarse *= HGCalParameters::k_ScaleToDDD;
0135 #ifdef EDM_ML_DEBUG
0136 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: Z " << zpos << ":" << zz << " R " << rfine << ":" << rcoarse;
0137 #endif
0138 return std::make_pair(rfine, rcoarse);
0139 }
0140
0141 double HGCalWaferType::areaPolygon(std::vector<double> const& x, std::vector<double> const& y) {
0142 double area = 0.0;
0143 int n = static_cast<int>(x.size());
0144 int j = n - 1;
0145 for (int i = 0; i < n; ++i) {
0146 area += ((x[j] + x[i]) * (y[i] - y[j]));
0147 j = i;
0148 }
0149 return (0.5 * area);
0150 }
0151
0152 std::pair<double, double> HGCalWaferType::intersection(
0153 int k1, int k2, std::vector<double> const& x, std::vector<double> const& y, double xpos, double ypos, double rr) {
0154 double slope = (x[k1] - x[k2]) / (y[k1] - y[k2]);
0155 double interc = x[k1] - slope * y[k1];
0156 double xx[2], yy[2], dist[2];
0157 double v1 = std::sqrt((slope * slope + 1) * rr * rr - (interc * interc));
0158 yy[0] = (-slope * interc + v1) / (1 + slope * slope);
0159 yy[1] = (-slope * interc - v1) / (1 + slope * slope);
0160 for (int i = 0; i < 2; ++i) {
0161 xx[i] = (slope * yy[i] + interc);
0162 dist[i] = ((xx[i] - xpos) * (xx[i] - xpos)) + ((yy[i] - ypos) * (yy[i] - ypos));
0163 }
0164 #ifdef EDM_ML_DEBUG
0165 edm::LogVerbatim("HGCalGeom") << "HGCalWaferType: InterSection " << dist[0] << ":" << xx[0] << ":" << yy[0] << " vs "
0166 << dist[1] << ":" << xx[1] << ":" << yy[1];
0167 #endif
0168 if (dist[0] > dist[1])
0169 return std::make_pair(xx[1], yy[1]);
0170 else
0171 return std::make_pair(xx[0], yy[0]);
0172 }
