CMS 3D CMS Logo
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 <string>
7 
8 //#define EDM_ML_DEBUG
9 HGCalGeomTools::HGCalGeomTools() : factor_(1.0 / std::sqrt(3.0)) {}
10 
11 void HGCalGeomTools::radius(double zf,
12  double zb,
13  std::vector<double> const& zFront1,
14  std::vector<double> const& rFront1,
15  std::vector<double> const& slope1,
16  std::vector<double> const& zFront2,
17  std::vector<double> const& rFront2,
18  std::vector<double> const& slope2,
19  int flag,
20  std::vector<double>& zz,
21  std::vector<double>& rin,
22  std::vector<double>& rout) {
23  double dz1(0), dz2(0);
24  auto zf1 = std::lower_bound(zFront1.begin(), zFront1.end(), zf);
25  if (zf1 != zFront1.begin())
26  --zf1;
27  if (((zf1 + 1) != zFront1.end()) && (std::abs(*(zf1 + 1) - zf) < tol_)) {
28  ++zf1;
29  dz1 = 2 * tol_;
30  }
31  auto zf2 = std::lower_bound(zFront2.begin(), zFront2.end(), zf);
32  if (zf2 != zFront2.begin())
33  --zf2;
34  auto zb1 = std::lower_bound(zFront1.begin(), zFront1.end(), zb);
35  if (zb1 != zFront1.begin())
36  --zb1;
37  if ((zb1 != zFront1.begin()) && (std::abs(*zb1 - zb) < tol_)) {
38  --zb1;
39  dz2 = -2 * tol_;
40  }
41  auto zb2 = std::lower_bound(zFront2.begin(), zFront2.end(), zb);
42  if (zb2 != zFront2.begin())
43  --zb2;
44 #ifdef EDM_ML_DEBUG
45  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius:zf " << zf << " : " << *zf1 << " : " << *zf2 << " zb " << zb
46  << " : " << *zb1 << " : " << *zb2 << " Flag " << flag;
47 #endif
48  if ((zf1 == zb1) && (zf2 == zb2)) {
49 #ifdef EDM_ML_DEBUG
50  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius:Try 1: " << zf << ":" << zb << " dz " << dz1 << ":" << dz2;
51 #endif
52  zz.emplace_back(zf);
53  rin.emplace_back(radius(zf + dz1, zFront1, rFront1, slope1));
54  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
55  zz.emplace_back(zb);
56  rin.emplace_back(radius(zb + dz2, zFront1, rFront1, slope1));
57  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
58  } else if (zf1 == zb1) {
59 #ifdef EDM_ML_DEBUG
60  double z1 = std::max(*zf2, *zb2);
61  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius:Try 2:" << zf << ":" << *zb2 << " (" << z1 << ") : " << zb;
62 #endif
63  zz.emplace_back(zf);
64  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
65  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
66  if (slope(*zb2, zFront2, slope2) < tol_) {
67  zz.emplace_back(*zb2);
68  rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
69  rout.emplace_back(radius(*zb2 - tol_, zFront2, rFront2, slope2));
70  }
71  if ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_)) {
72  zz.emplace_back(*zb2);
73  rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
74  rout.emplace_back(radius(*zb2, zFront2, rFront2, slope2));
75  }
76  zz.emplace_back(zb);
77  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
78  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
79  } else if (zf2 == zb2) {
80 #ifdef EDM_ML_DEBUG
81  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius:Try 3: " << zf << ":" << *zb1 << ":" << zb;
82 #endif
83  zz.emplace_back(zf);
84  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
85  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
86  if (slope(*zb1, zFront1, slope1) < tol_) {
87  zz.emplace_back(*zb1);
88  rin.emplace_back(radius(*zb1 - tol_, zFront1, rFront1, slope1));
89  rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
90  }
91  if ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_)) {
92  zz.emplace_back(*zb1);
93  rin.emplace_back(radius(*zb1, zFront1, rFront1, slope1));
94  rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
95  }
96  zz.emplace_back(zb);
97  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
98  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
99  } else {
100  double z1 = std::min(*zf2, *zb1);
101  double z2 = std::max(*zf2, *zb1);
102 #ifdef EDM_ML_DEBUG
103  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius:Try 4: " << zf << ":" << z1 << " : " << z2 << ":" << zb;
104 #endif
105  zz.emplace_back(zf);
106  rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
107  rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
108  zz.emplace_back(z1);
109  rin.emplace_back(radius(z1, zFront1, rFront1, slope1));
110  rout.emplace_back(radius(z1, zFront2, rFront2, slope2));
111  zz.emplace_back(z2);
112  rin.emplace_back(radius(z2, zFront1, rFront1, slope1));
113  rout.emplace_back(radius(z2, zFront2, rFront2, slope2));
114  zz.emplace_back(zb);
115  rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
116  rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
117  }
118  double rmin = *(std::min_element(rout.begin(), rout.end()));
119  if (flag > 1) {
120  for (auto& rr : rout)
121  rr = rmin;
122  }
123 #ifdef EDM_ML_DEBUG
124  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius has " << zz.size() << " sections: " << rmin;
125  for (unsigned int k = 0; k < zz.size(); ++k)
126  edm::LogVerbatim("HGCalGeom") << "[" << k << "] Z = " << zz[k] << " R = " << rin[k] << ":" << rout[k];
127 #endif
128 }
129 
130 double HGCalGeomTools::radius(double z,
131  std::vector<double> const& zFront,
132  std::vector<double> const& rFront,
133  std::vector<double> const& slope) {
134  auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
135  if (itrz != zFront.begin())
136  --itrz;
137  unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
138  if ((ik + 1) < zFront.size() && std::abs(z - zFront[ik + 1]) < tol_)
139  ++ik;
140  double r = rFront[ik] + (z - zFront[ik]) * slope[ik];
141 #ifdef EDM_ML_DEBUG
142  edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << " k " << ik << " R " << r;
143 #endif
144  return r;
145 }
146 
147 double HGCalGeomTools::radius(
148  double z, int layer0, int layerf, std::vector<double> const& zFront, std::vector<double> const& rFront) {
149  double r = rFront[0];
150 #ifdef EDM_ML_DEBUG
151  unsigned int ik(0);
152 #endif
153  for (unsigned int k = 0; k < rFront.size(); ++k) {
154  int k1 = layerf - layer0 + (int)(k);
155  if (k1 < (int)(zFront.size())) {
156  r = rFront[k];
157 #ifdef EDM_ML_DEBUG
158  ik = k;
159 #endif
160  if (z < zFront[k1] + tol_)
161  break;
162  }
163  }
164 #ifdef EDM_ML_DEBUG
165  edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << ":" << ik << " R " << r;
166 #endif
167  return r;
168 }
169 
170 std::pair<double, double> HGCalGeomTools::shiftXY(int waferPosition, double waferSize) const {
171  double dx(0), dy(0);
172  switch (waferPosition) {
173  case (HGCalTypes::CornerCenterYp): {
174  dy = factor_ * waferSize;
175  break;
176  }
177  case (HGCalTypes::CornerCenterYm): {
178  dy = -factor_ * waferSize;
179  break;
180  }
181  case (HGCalTypes::CornerCenterXp): {
182  dx = factor_ * waferSize;
183  break;
184  }
185  case (HGCalTypes::CornerCenterXm): {
186  dx = -factor_ * waferSize;
187  break;
188  }
189  }
190 #ifdef EDM_ML_DEBUG
191  edm::LogVerbatim("HGCalGeom") << "Shift for " << waferPosition << " is (" << dx << ":" << dy << ")";
192 #endif
193  return std::make_pair(dx, dy);
194 }
195 
196 double HGCalGeomTools::slope(double z, std::vector<double> const& zFront, std::vector<double> const& slope) {
197  auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
198  if (itrz != zFront.begin())
199  --itrz;
200  unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
201  // if (ik < zFront.size() && std::abs(z-zFront[ik+1]) < tol_) ++ik;
202 #ifdef EDM_ML_DEBUG
203  edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::slope:z " << z << " k " << ik << " Slope " << slope[ik];
204 #endif
205  return slope[ik];
206 }
207 
208 std::pair<double, double> HGCalGeomTools::zradius(double z1,
209  double z2,
210  std::vector<double> const& zF,
211  std::vector<double> const& rF) {
212  double z(-1), r(-1);
213  for (unsigned int k = 0; k < rF.size(); ++k) {
214  if ((z1 > zF[k] - tol_) && (z2 < zF[k] + tol_)) {
215  z = zF[k];
216  r = rF[k];
217  break;
218  }
219  }
220  return std::make_pair(z, r);
221 }
222 
223 std::pair<int32_t, int32_t> HGCalGeomTools::waferCorner(
224  double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug) {
225  double xc[HGCalParameters::k_CornerSize], yc[HGCalParameters::k_CornerSize];
226  xc[0] = xpos;
227  yc[0] = ypos + R;
228  xc[1] = xpos - r;
229  yc[1] = ypos + 0.5 * R;
230  if (oldBug) {
231  xc[2] = xpos + r;
232  yc[2] = ypos - 0.5 * R;
233  } else {
234  xc[2] = xpos - r;
235  yc[2] = ypos - 0.5 * R;
236  }
237  xc[3] = xpos;
238  yc[3] = ypos - R;
239  xc[4] = xpos + r;
240  yc[4] = ypos - 0.5 * R;
241  xc[5] = xpos + r;
242  yc[5] = ypos + 0.5 * R;
243  int32_t nCorner(0), firstCorner(-1), firstMiss(-1);
244 #ifdef EDM_ML_DEBUG
245  std::vector<uint32_t> corners;
246 #endif
247  for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
248  double rpos = sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
249  if ((rpos <= rMax) && (rpos >= rMin)) {
250 #ifdef EDM_ML_DEBUG
251  corners.emplace_back(k);
252 #endif
253  if (firstCorner < 0)
254  firstCorner = k;
255  ++nCorner;
256  } else {
257  if (firstMiss < 0)
258  firstMiss = k;
259  }
260  }
261  if ((nCorner > 1) && (firstCorner == 0) && (firstMiss < nCorner)) {
262  firstCorner = firstMiss + HGCalParameters::k_CornerSize - nCorner;
263  }
264 #ifdef EDM_ML_DEBUG
265  edm::LogVerbatim("HGCalGeom") << "waferCorner:: R " << rMin << ":" << rMax << nCorner << " corners; first corner "
266  << firstCorner;
267  for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
268  double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
269  std::string ok = (std::find(corners.begin(), corners.end(), k) != corners.end()) ? " In" : " Out";
270  edm::LogVerbatim("HGCalGeom") << "Corner[" << k << "] x " << xc[k] << " y " << yc[k] << " R " << rpos << ok;
271  }
272 #endif
273  return std::make_pair(nCorner, firstCorner);
274 }
photonAnalyzer_cfi.rMax
rMax
Definition: photonAnalyzer_cfi.py:91
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:223
geometryCSVtoXML.zz
zz
Definition: geometryCSVtoXML.py:19
MessageLogger.h
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:14
HGCalGeomTools::slope
static double slope(double z, std::vector< double > const &zFront, std::vector< double > const &slope)
Definition: HGCalGeomTools.cc:196
HGCalTypes::CornerCenterXm
Definition: HGCalTypes.h:44
min
T min(T a, T b)
Definition: MathUtil.h:58
findQualityFiles.rr
string rr
Definition: findQualityFiles.py:185
HGCalGeomTools::factor_
double factor_
Definition: HGCalGeomTools.h:47
HGCalParameters::k_CornerSize
static constexpr uint32_t k_CornerSize
Definition: HGCalParameters.h:38
HGCalTypes::CornerCenterYm
Definition: HGCalTypes.h:42
HGCalGeomTools::shiftXY
std::pair< double, double > shiftXY(int waferPosition, double waferSize) const
Definition: HGCalGeomTools.cc:170
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
testProducerWithPsetDescEmpty_cfi.z2
z2
Definition: testProducerWithPsetDescEmpty_cfi.py:41
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:208
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
DDAxes::z
dqmdumpme.k
k
Definition: dqmdumpme.py:60
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HGCalTypes::CornerCenterYp
Definition: HGCalTypes.h:41
pfDeepBoostedJetPreprocessParams_cfi.lower_bound
lower_bound
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:15
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
createfilelist.int
int
Definition: createfilelist.py:10
PVValHelper::dy
Definition: PVValidationHelpers.h:49
alignCSCRings.r
r
Definition: alignCSCRings.py:93
HGCalGeomTools.h
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:11
std
Definition: JetResolutionObject.h:76
HGCalGeomTools::HGCalGeomTools
HGCalGeomTools()
Definition: HGCalGeomTools.cc:9
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
HGCalTypes::CornerCenterXp
Definition: HGCalTypes.h:43
photonAnalyzer_cfi.rMin
rMin
Definition: photonAnalyzer_cfi.py:90
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
HGCalTypes.h
edm::Log
Definition: MessageLogger.h:70
dttmaxenums::R
Definition: DTTMax.h:29
HGCalParameters.h
HGCalGeomTools::tol_
static constexpr double tol_
Definition: HGCalGeomTools.h:46
PVValHelper::dx
Definition: PVValidationHelpers.h:48
RemoveAddSevLevel.flag
flag
Definition: RemoveAddSevLevel.py:116
Generated for CMSSW Reference Manual by doxygen 1.8.16