CMS 3D CMS Logo

HGCalWaferType.cc
Go to the documentation of this file.
4 
5 //#define EDM_ML_DEBUG
6 
7 HGCalWaferType::HGCalWaferType(const std::vector<double>& rad100,
8  const std::vector<double>& rad200,
9  double waferSize, double zMin, int choice,
10  unsigned int cornerCut, double cutArea)
11  : rad100_(rad100),
12  rad200_(rad200),
13  waferSize_(waferSize),
14  zMin_(zMin),
15  choice_(choice),
16  cutValue_(cornerCut),
17  cutFracArea_(cutArea) {
18  r_ = 0.5 * waferSize_;
19  R_ = sqrt3_ * waferSize_;
20 #ifdef EDM_ML_DEBUG
21  edm::LogVerbatim("HGCalGeom")
22  << "HGCalWaferType: initialized with waferR's " << waferSize_ << ":"
23  << r_ << ":" << R_ << " Choice " << choice_ << " Cuts " << cutValue_
24  << ":" << cutFracArea_ << " zMin " << zMin_ << " with " << rad100_.size()
25  << ":" << rad200_.size() << " parameters for R:";
26  for (unsigned k = 0; k < rad100_.size(); ++k)
27  edm::LogVerbatim("HGCalGeom")
28  << "[" << k << "] 100:200 " << rad100_[k] << " 200:300 " << rad200_[k];
29 #endif
30 }
31 
33 
34 int HGCalWaferType::getType(double xpos, double ypos, double zpos) {
35  std::vector<double> xc(HGCalParameters::k_CornerSize, 0);
36  std::vector<double> yc(HGCalParameters::k_CornerSize, 0);
37  xc[0] = xpos + r_;
38  yc[0] = ypos + 0.5 * R_;
39  xc[1] = xpos;
40  yc[1] = ypos + R_;
41  xc[2] = xpos - r_;
42  yc[2] = ypos + 0.5 * R_;
43  xc[3] = xpos - r_;
44  yc[3] = ypos - 0.5 * R_;
45  xc[4] = xpos;
46  yc[4] = ypos - R_;
47  xc[5] = xpos + r_;
48  yc[5] = ypos - 0.5 * R_;
49  std::pair<double, double> rv = rLimits(zpos);
50  std::vector<int> fine, coarse;
51  for (unsigned int k = 0; k < HGCalParameters::k_CornerSize; ++k) {
52  double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
53  if (rpos <= rv.first)
54  fine.emplace_back(k);
55  else if (rpos <= rv.second)
56  coarse.emplace_back(k);
57  }
58  int type(-2);
59  double fracArea(0);
60  if (choice_ == 1) {
61  if (fine.size() >= cutValue_)
62  type = 0;
63  else if (coarse.size() >= cutValue_)
64  type = 1;
65  else
66  type = 2;
67  } else {
68  if (fine.size() >= 4)
69  type = 0;
70  else if (coarse.size() >= 4 && fine.size() <= 1)
71  type = 1;
72  else if (coarse.size() < 2 && fine.empty())
73  type = 2;
74  else if (!fine.empty())
75  type = -1;
76  if (type < 0) {
77  unsigned int kmax = (type == -1) ? fine.size() : coarse.size();
78  std::vector<double> xcn, ycn;
79  for (unsigned int k = 0; k < kmax; ++k) {
80  unsigned int k1 = (type == -1) ? fine[k] : coarse[k];
81  unsigned int k2 = (k1 == xc.size() - 1) ? 0 : k1 + 1;
82  bool ok = ((type == -1)
83  ? (std::find(fine.begin(), fine.end(), k2) != fine.end())
84  : (std::find(coarse.begin(), coarse.end(), k2) !=
85  coarse.end()));
86  xcn.emplace_back(xc[k1]);
87  ycn.emplace_back(yc[k1]);
88  if (!ok) {
89  double rr = (type == -1) ? rv.first : rv.second;
90  std::pair<double, double> xy =
91  intersection(k1, k2, xc, yc, xpos, ypos, rr);
92  xcn.emplace_back(xy.first);
93  ycn.emplace_back(xy.second);
94  }
95  }
96  fracArea = areaPolygon(xcn, ycn) / areaPolygon(xc, yc);
97  type = (fracArea > cutFracArea_) ? -(1 + type) : -type;
98  }
99  }
100 #ifdef EDM_ML_DEBUG
101  edm::LogVerbatim("HGCalGeom")
102  << "HGCalWaferType: position " << xpos << ":" << ypos << ":" << zpos
103  << " R "
104  << ":" << rv.first << ":" << rv.second << " corners|type " << fine.size()
105  << ":" << coarse.size() << ":" << fracArea << ":" << type;
106 #endif
107  return type;
108 }
109 
110 std::pair<double, double> HGCalWaferType::rLimits(double zpos) {
111  double zz = std::abs(zpos);
112  if (zz < zMin_) zz = zMin_;
114  double rfine = rad100_[0];
115  double rcoarse = rad200_[0];
116  for (int i = 1; i < 5; ++i) {
117  rfine *= zz;
118  rfine += rad100_[i];
119  rcoarse *= zz;
120  rcoarse += rad200_[i];
121  }
122  return std::pair<double, double>(rfine * HGCalParameters::k_ScaleToDDD,
124 }
125 
126 double HGCalWaferType::areaPolygon(std::vector<double> const& x,
127  std::vector<double> const& y) {
128  double area = 0.0;
129  int n = x.size();
130  int j = n - 1;
131  for (int i = 0; i < n; ++i) {
132  area += ((x[j] + x[i]) * (y[i] - y[j]));
133  j = i;
134  }
135  return (0.5 * area);
136 }
137 
138 std::pair<double, double> HGCalWaferType::intersection(
139  int k1, int k2, std::vector<double> const& x, std::vector<double> const& y,
140  double xpos, double ypos, double rr) {
141  double slope = (x[k1] - x[k2]) / (y[k1] - y[k2]);
142  double interc = x[k1] - slope * y[k1];
143  double xx[2], yy[2], dist[2];
144  double v1 = std::sqrt((slope * slope + 1) * rr * rr - (interc * interc));
145  yy[0] = (-slope * interc + v1) / (1 + slope * slope);
146  yy[1] = (-slope * interc - v1) / (1 + slope * slope);
147  for (int i = 0; i < 2; ++i) {
148  xx[i] = (slope * yy[i] + interc);
149  dist[i] =
150  ((xx[i] - xpos) * (xx[i] - xpos)) + ((yy[i] - ypos) * (yy[i] - ypos));
151  }
152  if (dist[0] > dist[1])
153  return std::pair<double, double>(xx[1], yy[1]);
154  else
155  return std::pair<double, double>(xx[0], yy[0]);
156 }
const double waferSize_
type
Definition: HCALResponse.h:21
const std::vector< double > rad200_
const unsigned int cutValue_
const double cutFracArea_
static const double slope[3]
const double sqrt3_
std::pair< double, double > rLimits(double zpos)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
const int choice_
int getType(double xpos, double ypos, double zpos)
double areaPolygon(std::vector< double > const &, std::vector< double > const &)
T sqrt(T t)
Definition: SSEVec.h:18
std::pair< double, double > intersection(int, int, std::vector< double > const &, std::vector< double > const &, double xp, double yp, double rr)
static double k_ScaleFromDDD
const std::vector< double > rad100_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HGCalWaferType(const std::vector< double > &rad100, const std::vector< double > &rad200, double waferSize, double zMin, int choice, unsigned int cutValue, double cutFracArea)
const double zMin_
int k[5][pyjets_maxn]
static uint32_t k_CornerSize
static double k_ScaleToDDD