CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HGCalGeomTools.cc
Go to the documentation of this file.
1 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
2 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
3 #include "FWCore/MessageLogger/interface/MessageLogger.h"
4 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
5 
6 #include <sstream>
7 #include <string>
8 
9 //#define EDM_ML_DEBUG
10 HGCalGeomTools::HGCalGeomTools() : factor_(1.0 / std::sqrt(3.0)) {}
11 
12 void HGCalGeomTools::radius(double zf,
13  double zb,
14  std::vector<double> const& zFront1,
15  std::vector<double> const& rFront1,
16  std::vector<double> const& slope1,
17  std::vector<double> const& zFront2,
18  std::vector<double> const& rFront2,
19  std::vector<double> const& slope2,
20  int flag,
21  std::vector<double>& zz,
22  std::vector<double>& rin,
23  std::vector<double>& rout) {
24 #ifdef EDM_ML_DEBUG
25  std::ostringstream st1;
26  st1 << "Z = " << zf << ":" << zb << "; zFront1 =";
27  for (const auto& v : zFront1)
28  st1 << " " << v;
29  st1 << "; rFront1 =";
30  for (const auto& v : rFront1)
31  st1 << " " << v;
32  st1 << "; slope1 =";
33  for (const auto& v : slope1)
34  st1 << " " << v;
35  st1 << "; zFront2 =";
36  for (const auto& v : zFront2)
37  st1 << " " << v;
38  st1 << "; rFront2 =";
39  for (const auto& v : rFront2)
40  st1 << " " << v;
41  st1 << "; slope2 =";
42  for (const auto& v : slope2)
43  st1 << " " << v;
44  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX: " << st1.str();
45  double dz2f(0), dz2b(0);
46 #endif
47  double dz1f(0), dz1b(0);
48  auto zf1 = std::lower_bound(zFront1.begin(), zFront1.end(), zf);
49  if (zf1 != zFront1.begin())
50  --zf1;
51  if (((zf1 + 1) != zFront1.end()) && (std::abs(*(zf1 + 1) - zf) < tol_)) {
52  ++zf1;
53  dz1f = 2 * tol_;
54  }
55  auto zf2 = std::lower_bound(zFront2.begin(), zFront2.end(), zf);
56  if (zf2 != zFront2.begin())
57  --zf2;
58  if (((zf2 + 1) != zFront2.end()) && (std::abs(*(zf2 + 1) - zf) < tol_)) {
59  if (static_cast<int>(zf2 - zFront2.begin()) >= 2)
60  ++zf2;
61 #ifdef EDM_ML_DEBUG
62  dz2f = 2 * tol_;
63 #endif
64  }
65  auto zb1 = std::lower_bound(zFront1.begin(), zFront1.end(), zb);
66  if (zb1 != zFront1.begin())
67  --zb1;
68  if ((zb1 != zFront1.begin()) && (std::abs(*zb1 - zb) < tol_)) {
69  --zb1;
70  dz1b = -2 * tol_;
71  }
72  if (((zb1 + 1) != zFront1.end()) && (std::abs(*(zb1 + 1) - zb) < tol_)) {
73  dz1b = -2 * tol_;
74  }
75  auto zb2 = std::lower_bound(zFront2.begin(), zFront2.end(), zb);
76  if (zb2 != zFront2.begin()) {
77  --zb2;
78  }
79  if ((zb2 != zFront2.begin()) && (std::abs(*zb2 - zb) < tol_)) {
80  if (static_cast<int>(zb2 - zFront2.begin()) > 2)
81  --zb2;
82 #ifdef EDM_ML_DEBUG
83  dz2b = -2 * tol_;
84 #endif
85  }
86 #ifdef EDM_ML_DEBUG
87  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:zf " << zf << " : "
88  << static_cast<int>(zf1 - zFront1.begin()) << ":" << *zf1 << " : "
89  << static_cast<int>(zf2 - zFront2.begin()) << ":" << *zf2 << " zb " << zb << ":"
90  << static_cast<int>(zb1 - zFront1.begin()) << " : " << *zb1 << " : "
91  << static_cast<int>(zb2 - zFront2.begin()) << ":" << *zb2 << " Flag " << flag << ":"
92  << (zf1 == zb1) << ":" << (zf2 == zb2) << " dz " << dz1f << ":" << dz2f << ":" << dz1b
93  << ":" << dz2b;
94 #endif
95  if ((zf1 == zb1) && (zf2 == zb2)) {
96 #ifdef EDM_ML_DEBUG
97  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 1: " << zf << ":" << zb;
98 #endif
99  zz.emplace_back(zf);
100  rin.emplace_back(radius(zf + dz1f, zFront1, rFront1, slope1));
101  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
102  zz.emplace_back(zb);
103  rin.emplace_back(radius(zb + dz1b, zFront1, rFront1, slope1));
104  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
105  } else if (zf1 == zb1) {
106 #ifdef EDM_ML_DEBUG
107  double z1 = std::max(*zf2, *zb2);
108  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 2:" << zf << ":" << *zb2 << " (" << z1 << ") : " << zb
109  << " Test " << (slope(*zb2, zFront2, slope2) < tol_) << ":"
110  << ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_));
111 #endif
112  zz.emplace_back(zf);
113  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
114  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
115  if (slope(*zb2, zFront2, slope2) < tol_) {
116  zz.emplace_back(*zb2);
117  rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
118  rout.emplace_back(radius(*zb2 - tol_, zFront2, rFront2, slope2));
119  }
120  if ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_)) {
121  zz.emplace_back(*zb2);
122  rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
123  rout.emplace_back(radius(*zb2, zFront2, rFront2, slope2));
124  }
125  zz.emplace_back(zb);
126  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
127  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
128  } else if (zf2 == zb2) {
129 #ifdef EDM_ML_DEBUG
130  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 3: " << zf << ":" << *zb1 << ":" << zb << " Test "
131  << (slope(*zb1, zFront1, slope1) < tol_) << ":"
132  << ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_));
133 #endif
134  zz.emplace_back(zf);
135  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
136  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
137  if (slope(*zb1, zFront1, slope1) < tol_) {
138  zz.emplace_back(*zb1);
139  rin.emplace_back(radius(*zb1 - tol_, zFront1, rFront1, slope1));
140  rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
141  }
142  if ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_)) {
143  zz.emplace_back(*zb1);
144  rin.emplace_back(radius(*zb1, zFront1, rFront1, slope1));
145  rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
146  }
147  zz.emplace_back(zb);
148  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
149  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
150  } else {
151  double z1 = std::min(*zf2, *zb1);
152  double z2 = std::max(*zf2, *zb1);
153 #ifdef EDM_ML_DEBUG
154  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 4: " << zf << ":" << z1 << " : " << z2 << ":" << zb;
155 #endif
156  zz.emplace_back(zf);
157  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
158  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
159  zz.emplace_back(z1);
160  rin.emplace_back(radius(z1, zFront1, rFront1, slope1));
161  rout.emplace_back(radius(z1, zFront2, rFront2, slope2));
162  zz.emplace_back(z2);
163  rin.emplace_back(radius(z2, zFront1, rFront1, slope1));
164  rout.emplace_back(radius(z2, zFront2, rFront2, slope2));
165  zz.emplace_back(zb);
166  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
167  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
168  }
169  double rmin = *(std::min_element(rout.begin(), rout.end()));
170  if (flag > 1) {
171  for (auto& rr : rout)
172  rr = rmin;
173  }
174 #ifdef EDM_ML_DEBUG
175  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX has " << zz.size() << " sections: " << rmin;
176  for (unsigned int k = 0; k < zz.size(); ++k)
177  edm::LogVerbatim("HGCalGeom") << "[" << k << "] Z = " << zz[k] << " R = " << rin[k] << ":" << rout[k];
178 #endif
179 }
180 
181 double HGCalGeomTools::radius(double z,
182  std::vector<double> const& zFront,
183  std::vector<double> const& rFront,
184  std::vector<double> const& slope) {
185 #ifdef EDM_ML_DEBUG
186  std::ostringstream st1;
187  st1 << "Z = " << z << "; zFront =";
188  for (const auto& v : zFront)
189  st1 << " " << v;
190  st1 << "; rFront =";
191  for (const auto& v : rFront)
192  st1 << " " << v;
193  st1 << "; slope =";
194  for (const auto& v : slope)
195  st1 << " " << v;
196  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius: " << st1.str();
197 #endif
198  auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
199  if (itrz != zFront.begin())
200  --itrz;
201  unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
202  if ((ik + 1) < zFront.size() && std::abs(z - zFront[ik + 1]) < tol_)
203  ++ik;
204  double r = rFront[ik] + (z - zFront[ik]) * slope[ik];
205 #ifdef EDM_ML_DEBUG
206  edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << " k " << ik << " R " << r;
207 #endif
208  return r;
209 }
210 
211 double HGCalGeomTools::radius(
212  double z, int layer0, int layerf, std::vector<double> const& zFront, std::vector<double> const& rFront) {
213  double r = rFront[0];
214 #ifdef EDM_ML_DEBUG
215  unsigned int ik(0);
216 #endif
217  for (unsigned int k = 0; k < rFront.size(); ++k) {
218  int k1 = layerf - layer0 + (int)(k);
219  if (k1 < (int)(zFront.size())) {
220  r = rFront[k];
221 #ifdef EDM_ML_DEBUG
222  ik = k;
223 #endif
224  if (z < zFront[k1] + tol_)
225  break;
226  }
227  }
228 #ifdef EDM_ML_DEBUG
229  edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << ":" << ik << " R " << r;
230 #endif
231  return r;
232 }
233 
234 std::pair<double, double> HGCalGeomTools::shiftXY(int waferPosition, double waferSize) const {
235  double dx(0), dy(0);
236  switch (waferPosition) {
237  case (HGCalTypes::CornerCenterYp): {
238  dy = factor_ * waferSize;
239  break;
240  }
241  case (HGCalTypes::CornerCenterYm): {
242  dy = -factor_ * waferSize;
243  break;
244  }
245  case (HGCalTypes::CornerCenterXp): {
246  dx = factor_ * waferSize;
247  break;
248  }
249  case (HGCalTypes::CornerCenterXm): {
250  dx = -factor_ * waferSize;
251  break;
252  }
253  }
254 #ifdef EDM_ML_DEBUG
255  edm::LogVerbatim("HGCalGeom") << "Shift for " << waferPosition << " is (" << dx << ":" << dy << ")";
256 #endif
257  return std::make_pair(dx, dy);
258 }
259 
260 double HGCalGeomTools::slope(double z, std::vector<double> const& zFront, std::vector<double> const& slope) {
261  auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
262  if (itrz != zFront.begin())
263  --itrz;
264  unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
265  // if (ik < zFront.size() && std::abs(z-zFront[ik+1]) < tol_) ++ik;
266 #ifdef EDM_ML_DEBUG
267  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::slope:z " << z << " k " << ik << " Slope " << slope[ik];
268 #endif
269  return slope[ik];
270 }
271 
272 std::pair<double, double> HGCalGeomTools::zradius(double z1,
273  double z2,
274  std::vector<double> const& zF,
275  std::vector<double> const& rF) {
276  double z(-1), r(-1);
277  for (unsigned int k = 0; k < rF.size(); ++k) {
278  if ((z1 > zF[k] - tol_) && (z2 < zF[k] + tol_)) {
279  z = zF[k];
280  r = rF[k];
281  break;
282  }
283  }
284  return std::make_pair(z, r);
285 }
286 
287 std::pair<int32_t, int32_t> HGCalGeomTools::waferCorner(
288  double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug) {
289  double xc[HGCalParameters::k_CornerSize], yc[HGCalParameters::k_CornerSize];
290  xc[0] = xpos;
291  yc[0] = ypos + R;
292  xc[1] = xpos - r;
293  yc[1] = ypos + 0.5 * R;
294  if (oldBug) {
295  xc[2] = xpos + r;
296  yc[2] = ypos - 0.5 * R;
297  } else {
298  xc[2] = xpos - r;
299  yc[2] = ypos - 0.5 * R;
300  }
301  xc[3] = xpos;
302  yc[3] = ypos - R;
303  xc[4] = xpos + r;
304  yc[4] = ypos - 0.5 * R;
305  xc[5] = xpos + r;
306  yc[5] = ypos + 0.5 * R;
307  int32_t nCorner(0), firstCorner(-1), firstMiss(-1);
308 #ifdef EDM_ML_DEBUG
309  std::vector<uint32_t> corners;
310 #endif
311  for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
312  double rpos = sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
313  if ((rpos <= rMax) && (rpos >= rMin)) {
314 #ifdef EDM_ML_DEBUG
315  corners.emplace_back(k);
316 #endif
317  if (firstCorner < 0)
318  firstCorner = k;
319  ++nCorner;
320  } else {
321  if (firstMiss < 0)
322  firstMiss = k;
323  }
324  }
325  if ((nCorner > 1) && (firstCorner == 0) && (firstMiss < nCorner)) {
326  firstCorner = firstMiss + HGCalParameters::k_CornerSize - nCorner;
327  }
328 #ifdef EDM_ML_DEBUG
329  edm::LogVerbatim("HGCalGeom") << "waferCorner:: R " << rMin << ":" << rMax << nCorner << " corners; first corner "
330  << firstCorner;
331  for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
332  double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
333  std::string ok = (std::find(corners.begin(), corners.end(), k) != corners.end()) ? " In" : " Out";
334  edm::LogVerbatim("HGCalGeom") << "Corner[" << k << "] x " << xc[k] << " y " << yc[k] << " R " << rpos << ok;
335  }
336 #endif
337  return std::make_pair(nCorner, firstCorner);
338 }
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
HGCalGeomTools::radius
static void radius(double zf, double zb, std::vector< double > const &zFront1, std::vector< double > const &rFront1, std::vector< double > const &slope1, std::vector< double > const &zFront2, std::vector< double > const &rFront2, std::vector< double > const &slope2, int flag, std::vector< double > &zz, std::vector< double > &rin, std::vector< double > &rout)
Definition: HGCalGeomTools.cc:12
HGCalGeomTools::tol_
static constexpr double tol_
Definition: HGCalGeomTools.h:47
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
HGCalGeomTools::slope
static double slope(double z, std::vector< double > const &zFront, std::vector< double > const &slope)
Definition: HGCalGeomTools.cc:260
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:14
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HGCalParameters::k_CornerSize
static constexpr uint32_t k_CornerSize
Definition: HGCalParameters.h:44
HGCalGeomTools::waferCorner
static std::pair< int32_t, int32_t > waferCorner(double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug=false)
Definition: HGCalGeomTools.cc:287
HGCalGeomTools::zradius
static std::pair< double, double > zradius(double z1, double z2, std::vector< double > const &zFront, std::vector< double > const &rFront)
Definition: HGCalGeomTools.cc:272
HGCalGeomTools::shiftXY
std::pair< double, double > shiftXY(int waferPosition, double waferSize) const
Definition: HGCalGeomTools.cc:234
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
cuda_std::lower_bound
__host__ __device__ constexpr RandomIt lower_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
Definition: cudastdAlgorithm.h:27
Generated for CMSSW Reference Manual by doxygen 1.8.5