Back to home page

Project CMSSW displayed by LXR

 
 

    


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