Back to home page

Project CMSSW displayed by LXR

 
 

    


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