Belle II Software  release-05-02-19
GeometryData.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Kirill Chilikin *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <klm/eklm/geometry/GeometryData.h>
13 
14 /* KLM headers. */
15 #include <klm/eklm/geometry/Circle2D.h>
16 #include <klm/eklm/geometry/Line2D.h>
17 
18 /* Belle 2 headers. */
19 #include <framework/database/DBObjPtr.h>
20 #include <framework/database/Database.h>
21 #include <framework/logging/Logger.h>
22 
23 /* CLHEP headers. */
24 #include <CLHEP/Geometry/Point3D.h>
25 #include <CLHEP/Units/SystemOfUnits.h>
26 
27 /* C++ headers. */
28 #include <string>
29 
30 using namespace Belle2;
31 
32 static const char c_MemErr[] = "Memory allocation error.";
33 
34 const EKLM::GeometryData&
35 EKLM::GeometryData::Instance(enum DataSource dataSource, const GearDir* gearDir)
36 {
37  static EKLM::GeometryData gd(dataSource, gearDir);
38  return gd;
39 }
40 
46 static void readPositionData(EKLMGeometry::ElementPosition* epos, GearDir* gd)
47 {
48  epos->setX(gd->getLength("X") * CLHEP::cm);
49  epos->setY(gd->getLength("Y") * CLHEP::cm);
50  epos->setZ(gd->getLength("Z") * CLHEP::cm);
51 }
52 
58 static void readSizeData(EKLMGeometry::ElementPosition* epos, GearDir* gd)
59 {
60  epos->setInnerR(gd->getLength("InnerR") * CLHEP::cm);
61  epos->setOuterR(gd->getLength("OuterR") * CLHEP::cm);
62  epos->setLength(gd->getLength("Length") * CLHEP::cm);
63 }
64 
70 static void readSectorSupportGeometry(
71  EKLMGeometry::SectorSupportGeometry* ssg, GearDir* gd)
72 {
73  ssg->setThickness(gd->getLength("Thickness") * CLHEP::cm);
74  ssg->setDeltaLY(gd->getLength("DeltaLY") * CLHEP::cm);
75  ssg->setCornerX(gd->getLength("CornerX") * CLHEP::cm);
76  ssg->setCorner1LX(gd->getLength("Corner1LX") * CLHEP::cm);
77  ssg->setCorner1Width(gd->getLength("Corner1Width") * CLHEP::cm);
78  ssg->setCorner1Thickness(gd->getLength("Corner1Thickness") * CLHEP::cm);
79  ssg->setCorner1Z(gd->getLength("Corner1Z") * CLHEP::cm);
80  ssg->setCorner2LX(gd->getLength("Corner2LX") * CLHEP::cm);
81  ssg->setCorner2LY(gd->getLength("Corner2LY") * CLHEP::cm);
82  ssg->setCorner2Thickness(gd->getLength("Corner2Thickness") * CLHEP::cm);
83  ssg->setCorner2Z(gd->getLength("Corner2Z") * CLHEP::cm);
84  ssg->setCorner3LX(gd->getLength("Corner3LX") * CLHEP::cm);
85  ssg->setCorner3LY(gd->getLength("Corner3LY") * CLHEP::cm);
86  ssg->setCorner3Thickness(gd->getLength("Corner3Thickness") * CLHEP::cm);
87  ssg->setCorner3Z(gd->getLength("Corner3Z") * CLHEP::cm);
88  ssg->setCorner4LX(gd->getLength("Corner4LX") * CLHEP::cm);
89  ssg->setCorner4LY(gd->getLength("Corner4LY") * CLHEP::cm);
90  ssg->setCorner4Thickness(gd->getLength("Corner4Thickness") * CLHEP::cm);
91  ssg->setCorner4Z(gd->getLength("Corner4Z") * CLHEP::cm);
92 }
93 
99 static void readShieldDetailGeometry(
100  EKLMGeometry::ShieldDetailGeometry* sdg, GearDir* gd)
101 {
102  int i, n;
103  std::string name;
104  EKLMGeometry::Point p;
105  sdg->setLengthX(gd->getLength("LengthX") * CLHEP::cm);
106  sdg->setLengthY(gd->getLength("LengthY") * CLHEP::cm);
107  n = gd->getNumberNodes("Point");
108  sdg->setNPoints(n);
109  for (i = 0; i < n; i++) {
110  GearDir point(*gd);
111  name = "/Point[" + std::to_string(i + 1) + "]";
112  point.append(name.c_str());
113  p.setX(point.getLength("X") * CLHEP::cm);
114  p.setY(point.getLength("Y") * CLHEP::cm);
115  sdg->setPoint(i, p);
116  }
117 }
118 
119 void EKLM::GeometryData::calculateSectorSupportGeometry()
120 {
121  HepGeom::Point3D<double> p;
122  Line2D line23Outer(0, m_SectorSupportPosition.getY(), 1, 0);
123  Line2D line23Inner(0, m_SectorSupportPosition.getY() +
124  m_SectorSupportGeometry.getThickness(), 1, 0);
125  Line2D line23Prism(0, m_SectorSupportPosition.getY() +
126  m_SectorSupportGeometry.getThickness() +
127  m_SectorSupportGeometry.getCorner3LY(), 1, 0);
128  Line2D line41Outer(m_SectorSupportPosition.getX(), 0, 0, 1);
129  Line2D line41Inner(m_SectorSupportPosition.getX() +
130  m_SectorSupportGeometry.getThickness(), 0, 0, 1);
131  Line2D line41Prism(m_SectorSupportPosition.getX() +
132  m_SectorSupportGeometry.getThickness() +
133  m_SectorSupportGeometry.getCorner4LX(), 0, 0, 1);
134  Line2D line41Corner1B(m_SectorSupportPosition.getX() +
135  m_SectorSupportGeometry.getCornerX(), 0, 0, 1);
136  Circle2D circleInnerOuter(0, 0, m_SectorSupportPosition.getInnerR());
137  Circle2D circleInnerInner(0, 0, m_SectorSupportPosition.getInnerR() +
138  m_SectorSupportGeometry.getThickness());
139  Circle2D circleOuterInner(0, 0, m_SectorSupportPosition.getOuterR() -
140  m_SectorSupportGeometry.getThickness());
141  Circle2D circleOuterOuter(0, 0, m_SectorSupportPosition.getOuterR());
142  HepGeom::Point3D<double> intersections[2];
143  /* Corner 1. */
144  p.setX(m_SectorSupportPosition.getX());
145  p.setY(m_SectorSupportPosition.getOuterR() -
146  m_SectorSupportGeometry.getDeltaLY());
147  p.setZ(0);
148  m_SectorSupportGeometry.setCorner1A(p);
149  line41Corner1B.findIntersection(circleOuterOuter, intersections);
150  m_SectorSupportGeometry.setCorner1B(intersections[1]);
151  m_SectorSupportGeometry.setCornerAngle(
152  atan2(m_SectorSupportGeometry.getCorner1B().y() -
153  m_SectorSupportGeometry.getCorner1A().y(),
154  m_SectorSupportGeometry.getCorner1B().x() -
155  m_SectorSupportGeometry.getCorner1A().x()) * CLHEP::rad);
156  p.setX(m_SectorSupportPosition.getX() +
157  m_SectorSupportGeometry.getThickness());
158  p.setY(m_SectorSupportGeometry.getCorner1A().y() -
159  m_SectorSupportGeometry.getThickness() *
160  (1.0 / cos(m_SectorSupportGeometry.getCornerAngle()) -
161  tan(m_SectorSupportGeometry.getCornerAngle())));
162  p.setZ(0);
163  m_SectorSupportGeometry.setCorner1AInner(p);
164  Line2D lineCorner1(m_SectorSupportGeometry.getCorner1AInner().x(),
165  m_SectorSupportGeometry.getCorner1AInner().y(),
166  m_SectorSupportGeometry.getCorner1B().x() -
167  m_SectorSupportGeometry.getCorner1A().x(),
168  m_SectorSupportGeometry.getCorner1B().y() -
169  m_SectorSupportGeometry.getCorner1A().y());
170  lineCorner1.findIntersection(circleOuterInner, intersections);
171  m_SectorSupportGeometry.setCorner1BInner(intersections[1]);
172  /* Corner 2. */
173  line23Inner.findIntersection(circleOuterInner, intersections);
174  m_SectorSupportGeometry.setCorner2Inner(intersections[1]);
175  /* Corner 3. */
176  line23Outer.findIntersection(circleInnerOuter, intersections);
177  m_SectorSupportGeometry.setCorner3(intersections[1]);
178  line23Inner.findIntersection(circleInnerInner, intersections);
179  m_SectorSupportGeometry.setCorner3Inner(intersections[1]);
180  line23Prism.findIntersection(circleInnerInner, intersections);
181  p.setX(intersections[1].x());
182  p.setY(m_SectorSupportPosition.getY() +
183  m_SectorSupportGeometry.getThickness());
184  p.setZ(0);
185  m_SectorSupportGeometry.setCorner3Prism(p);
186  /* Corner 4. */
187  line41Outer.findIntersection(circleInnerOuter, intersections);
188  m_SectorSupportGeometry.setCorner4(intersections[1]);
189  line41Inner.findIntersection(circleInnerInner, intersections);
190  m_SectorSupportGeometry.setCorner4Inner(intersections[1]);
191  line41Prism.findIntersection(circleInnerInner, intersections);
192  p.setX(m_SectorSupportPosition.getX() +
193  m_SectorSupportGeometry.getThickness());
194  p.setY(intersections[1].y());
195  p.setZ(0);
196  m_SectorSupportGeometry.setCorner4Prism(p);
197 }
198 
199 static bool compareLength(double a, double b)
200 {
201  return a < b;
202 }
203 
204 void EKLM::GeometryData::fillStripIndexArrays()
205 {
206  const char err[] = "Strip sorting algorithm error.";
207  int i;
208  double l;
209  std::vector<double> strips;
210  std::vector<double>::iterator it;
211  std::map<double, int> mapLengthStrip;
212  std::map<double, int> mapLengthStrip2;
213  std::map<double, int>::iterator itm;
214  for (i = 0; i < m_NStrips; i++) {
215  strips.push_back(m_StripPosition[i].getLength());
216  mapLengthStrip.insert(
217  std::pair<double, int>(m_StripPosition[i].getLength(), i));
218  }
219  sort(strips.begin(), strips.end(), compareLength);
220  l = strips[0];
221  m_nStripDifferent = 1;
222  for (it = strips.begin(); it != strips.end(); ++it) {
223  if ((*it) != l) {
224  l = (*it);
225  m_nStripDifferent++;
226  }
227  }
228  m_StripLenToAll = (int*)malloc(m_nStripDifferent * sizeof(int));
229  if (m_StripLenToAll == nullptr)
230  B2FATAL(c_MemErr);
231  i = 0;
232  l = strips[0];
233  itm = mapLengthStrip.find(l);
234  if (itm == mapLengthStrip.end())
235  B2FATAL(err);
236  m_StripLenToAll[i] = itm->second;
237  mapLengthStrip2.insert(std::pair<double, int>(l, i));
238  for (it = strips.begin(); it != strips.end(); ++it) {
239  if ((*it) != l) {
240  l = (*it);
241  i++;
242  itm = mapLengthStrip.find(l);
243  if (itm == mapLengthStrip.end())
244  B2FATAL(err);
245  m_StripLenToAll[i] = itm->second;
246  mapLengthStrip2.insert(std::pair<double, int>(l, i));
247  }
248  }
249  m_StripAllToLen = (int*)malloc(m_NStrips * sizeof(int));
250  if (m_StripAllToLen == nullptr)
251  B2FATAL(c_MemErr);
252  for (i = 0; i < m_NStrips; i++) {
253  itm = mapLengthStrip2.find(m_StripPosition[i].getLength());
254  if (itm == mapLengthStrip2.end())
255  B2FATAL(err);
256  m_StripAllToLen[i] = itm->second;
257  }
258 }
259 
260 void EKLM::GeometryData::readXMLDataStrips(const GearDir& gd)
261 {
262  int i;
263  std::string name;
264  GearDir Strips(gd);
265  Strips.append("/Strip");
266  m_StripGeometry.setWidth(Strips.getLength("Width") * CLHEP::cm);
267  m_StripGeometry.setThickness(Strips.getLength("Thickness") * CLHEP::cm);
268  m_StripGeometry.setGrooveDepth(Strips.getLength("GrooveDepth") * CLHEP::cm);
269  m_StripGeometry.setGrooveWidth(Strips.getLength("GrooveWidth") * CLHEP::cm);
270  m_StripGeometry.setNoScintillationThickness(
271  Strips.getLength("NoScintillationThickness") * CLHEP::cm);
272  m_StripGeometry.setRSSSize(Strips.getLength("RSSSize") * CLHEP::cm);
273  try {
274  m_StripPosition = new EKLMGeometry::ElementPosition[m_NStrips];
275  } catch (std::bad_alloc& ba) {
276  B2FATAL(c_MemErr);
277  }
278  for (i = 0; i < m_NStrips; i++) {
279  GearDir StripContent(Strips);
280  name = "/Strip[" + std::to_string(i + 1) + "]";
281  StripContent.append(name.c_str());
282  m_StripPosition[i].setLength(StripContent.getLength("Length") * CLHEP::cm);
283  m_StripPosition[i].setX(StripContent.getLength("X") * CLHEP::cm);
284  m_StripPosition[i].setY(StripContent.getLength("Y") * CLHEP::cm);
285  m_StripPosition[i].setZ(StripContent.getLength("Z") * CLHEP::cm);
286  }
287 }
288 
299 static void getDetailDxDy(HepGeom::Point3D<double>* points, int nPoints,
300  double r, double kx, double ky,
301  double& dx, double& dy)
302 {
303  int i;
304  /* Variable maxt is initialized to avoid a false-positive warning. */
305  /* cppcheck-suppress variableScope */
306  double a, b, c, d, t, maxt = 0, x1, y1, x2, y2, u;
307  bool intersection;
308  /*
309  * Contact by one of the detail points.
310  * Solve equation (x1 + kx * t)^2 + (y1 + ky * t)^2 = R^2,
311  * (kx^2 + ky^2) * t^2 + 2 * (kx * x1 + ky * y1) * t + x1^2 + y1^2 - r^2 = 0.
312  */
313  a = kx * kx + ky * ky;
314  intersection = false;
315  for (i = 0; i < nPoints; i++) {
316  x1 = points[i].x();
317  y1 = points[i].y();
318  b = 2.0 * (kx * x1 + ky * y1);
319  c = x1 * x1 + y1 * y1 - r * r;
320  d = b * b - 4.0 * a * c;
321  if (d >= 0) {
322  t = (-b + sqrt(d)) / (2.0 * a);
323  if (!intersection) {
324  intersection = true;
325  maxt = t;
326  } else {
327  if (t > maxt)
328  maxt = t;
329  }
330  }
331  }
332  if (!intersection)
333  B2FATAL("Shield layer geometry calculation error.");
334  /*
335  * Contact by one of the detail lines.
336  * Find t such as the equation
337  * (x1 + kx * t + (x2 - x1) * u)^2 + (y1 + ky * t + (y2 - y1) * u) = r^2
338  * has one solution relatively to u. Equation on t is
339  * 0 = ((x2 - x1) * (x1 + kx * t) + (y2 - y1) * (y1 + ky * t))^2 +
340  * ((x2 - x1)^2 + (y2 - y1)^2) * ((x1 + kx * t)^2 + (y1 + ky * t)^2 - r^2),
341  * t = (x1 * y2 - x2 * y1 +- r * sqrt((x2 - x1)^2 + (y2 - y1)^2)) /
342  * ((x2 - x1) * ky - (y2 - y1) * kx).
343  */
344  for (i = 0; i < nPoints; i++) {
345  x1 = points[i].x();
346  y1 = points[i].y();
347  if (i < nPoints - 1) {
348  x2 = points[i + 1].x();
349  y2 = points[i + 1].y();
350  } else {
351  x2 = points[0].x();
352  y2 = points[0].y();
353  }
354  a = (x2 - x1) * ky - (y2 - y1) * kx;
355  if (a == 0)
356  continue;
357  b = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
358  t = (x1 * y2 - x2 * y1 + r * sqrt(b)) / a;
359  /*
360  * Check whether intersection occurs between the translated points
361  * (x1 + kx * t, y1 + ky * t) and (x2 + kx * t, y2 + ky * t).
362  * (find solition of the orginal equation relatively to u for that).
363  */
364  u = -((x2 - x1) * (x1 + kx * t) + (y2 - y1) * (y2 + ky * t)) / b;
365  if (u < 0 || u > 1)
366  continue;
367  if (t > maxt)
368  maxt = t;
369  }
370  dx = kx * maxt;
371  dy = ky * maxt;
372 }
373 
379 static void EKLMPointToCLHEP(const EKLMGeometry::Point* pointEKLM,
380  HepGeom::Point3D<double>& pointCLHEP)
381 {
382  pointCLHEP.setX(pointEKLM->getX());
383  pointCLHEP.setY(pointEKLM->getY());
384  pointCLHEP.setZ(0);
385 }
386 
387 void EKLM::GeometryData::calculateShieldGeometry()
388 {
389  int i;
390  double r, l, dx, dy, xCenter, yCenter;
391  const double asqrt2 = 1.0 / sqrt(2.0);
392  HepGeom::Point3D<double> points[8];
393  const ShieldDetailGeometry* detailA = m_ShieldGeometry.getDetailA();
394  const ShieldDetailGeometry* detailB = m_ShieldGeometry.getDetailB();
395  const ShieldDetailGeometry* detailC = m_ShieldGeometry.getDetailC();
396  const ShieldDetailGeometry* detailD = m_ShieldGeometry.getDetailD();
397  r = m_SectorSupportPosition.getInnerR() +
398  m_SectorSupportGeometry.getThickness();
399  /* Detail A. */
400  EKLMPointToCLHEP(detailA->getPoint(0), points[0]);
401  EKLMPointToCLHEP(detailA->getPoint(1), points[1]);
402  EKLMPointToCLHEP(detailA->getPoint(2), points[2]);
403  EKLMPointToCLHEP(detailA->getPoint(3), points[3]);
404  EKLMPointToCLHEP(detailA->getPoint(4), points[4]);
405  points[5].setX(detailA->getLengthX());
406  points[5].setY(detailA->getLengthY());
407  points[5].setZ(0);
408  EKLMPointToCLHEP(detailA->getPoint(5), points[6]);
409  EKLMPointToCLHEP(detailA->getPoint(6), points[7]);
410  l = 0.5 * (detailA->getLengthX() + detailB->getLengthX());
411  xCenter = -asqrt2 * l;
412  yCenter = asqrt2 * l;
413  for (i = 0; i < 8; i++)
414  points[i] = HepGeom::Translate3D(xCenter, yCenter, 0) *
415  HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
416  HepGeom::Translate3D(-detailA->getLengthX() / 2,
417  -detailA->getLengthY() / 2, 0) *
418  points[i];
419  getDetailDxDy(points, 8, r, 1, 1, dx, dy);
420  m_ShieldGeometry.setDetailACenter(xCenter + dx, yCenter + dy);
421  /* Details B, D, E. */
422  points[0].setX(0);
423  points[0].setY(-detailD->getLengthY());
424  points[0].setZ(0);
425  points[1].setX(detailD->getLengthX());
426  points[1].setY(0);
427  points[1].setZ(0);
428  points[2].setX(detailB->getLengthX() - detailD->getLengthX());
429  points[2].setY(0);
430  points[2].setZ(0);
431  points[3].setX(detailB->getLengthX());
432  points[3].setY(-detailD->getLengthY());
433  points[3].setZ(0);
434  EKLMPointToCLHEP(detailB->getPoint(0), points[4]);
435  EKLMPointToCLHEP(detailB->getPoint(1), points[5]);
436  EKLMPointToCLHEP(detailB->getPoint(2), points[6]);
437  EKLMPointToCLHEP(detailB->getPoint(3), points[7]);
438  /* Detail B center coordinates before its shift are (0, 0). */
439  for (i = 0; i < 8; i++)
440  points[i] = HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
441  HepGeom::Translate3D(-detailB->getLengthX() / 2,
442  -detailB->getLengthY() / 2, 0) *
443  points[i];
444  getDetailDxDy(points, 8, r, 1, 1, dx, dy);
445  m_ShieldGeometry.setDetailBCenter(dx, dy);
446  /* Detail C. */
447  EKLMPointToCLHEP(detailC->getPoint(0), points[0]);
448  EKLMPointToCLHEP(detailC->getPoint(1), points[1]);
449  EKLMPointToCLHEP(detailC->getPoint(2), points[2]);
450  EKLMPointToCLHEP(detailC->getPoint(3), points[3]);
451  EKLMPointToCLHEP(detailC->getPoint(4), points[4]);
452  points[5].setX(detailC->getLengthX());
453  points[5].setY(detailC->getLengthY());
454  points[5].setZ(0);
455  EKLMPointToCLHEP(detailC->getPoint(5), points[6]);
456  EKLMPointToCLHEP(detailC->getPoint(6), points[7]);
457  l = 0.5 * (detailB->getLengthX() + detailC->getLengthX());
458  xCenter = asqrt2 * l;
459  yCenter = -asqrt2 * l;
460  for (i = 0; i < 8; i++)
461  points[i] = HepGeom::Translate3D(xCenter, yCenter, 0) *
462  HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
463  HepGeom::RotateY3D(180.0 * CLHEP::deg) *
464  HepGeom::Translate3D(-detailC->getLengthX() / 2,
465  -detailC->getLengthY() / 2, 0) *
466  points[i];
467  getDetailDxDy(points, 8, r, 1, 1, dx, dy);
468  m_ShieldGeometry.setDetailCCenter(xCenter + dx, yCenter + dy);
469 }
470 
471 void EKLM::GeometryData::readEndcapStructureGeometry(const GearDir& gd)
472 {
473  GearDir d(gd);
474  d.append("/EndcapStructure");
475  m_EndcapStructureGeometry.setPhi(d.getAngle("Phi") * CLHEP::rad);
476  m_EndcapStructureGeometry.setNSides(d.getInt("NSides"));
477 }
478 
479 void EKLM::GeometryData::initializeFromGearbox(const GearDir* gearDir)
480 {
481  int i, j, k;
482  std::string name;
483  ShieldDetailGeometry shieldDetailGeometry;
484  GearDir gd(*gearDir);
485  gd.append("/EKLM");
486  /* Numbers of elements. */
487  m_NSections = gd.getInt("NSections");
488  m_ElementNumbers->checkSection(m_NSections);
489  m_NLayers = gd.getInt("NLayers");
490  m_ElementNumbers->checkLayer(m_NLayers);
491  m_NDetectorLayers = new int[m_NSections];
492  m_NDetectorLayers[0] = gd.getInt("NDetectorLayersBackward");
493  checkDetectorLayerNumber(1, m_NDetectorLayers[0]);
494  if (m_NSections == 2) {
495  m_NDetectorLayers[1] = gd.getInt("NDetectorLayersForward");
496  checkDetectorLayerNumber(2, m_NDetectorLayers[1]);
497  }
498  m_NSectors = gd.getInt("NSectors");
499  m_ElementNumbers->checkSector(m_NSectors);
500  m_NPlanes = gd.getInt("NPlanes");
501  m_ElementNumbers->checkPlane(m_NPlanes);
502  m_NSegments = gd.getInt("NSegments");
503  m_ElementNumbers->checkSegment(m_NSegments);
504  m_NSegmentSupportElementsSector = (m_NSegments + 1) * m_NPlanes;
505  m_NStrips = gd.getInt("NStrips");
506  m_ElementNumbers->checkStrip(m_NStrips);
507  /* Geometry parameters. */
508  m_SolenoidZ = gd.getLength("SolenoidZ") * CLHEP::cm;
509  readEndcapStructureGeometry(gd);
510  GearDir section(gd);
511  section.append("/Section");
512  readPositionData(&m_SectionPosition, &section);
513  readSizeData(&m_SectionPosition, &section);
514  GearDir layer(gd);
515  layer.append("/Layer");
516  readSizeData(&m_LayerPosition, &layer);
517  m_LayerShiftZ = layer.getLength("ShiftZ") * CLHEP::cm;
518  GearDir sector(gd);
519  sector.append("/Sector");
520  readSizeData(&m_SectorPosition, &sector);
521  GearDir sectorSupport(gd);
522  sectorSupport.append("/SectorSupport");
523  readPositionData(&m_SectorSupportPosition, &sectorSupport);
524  readSizeData(&m_SectorSupportPosition, &sectorSupport);
525  readSectorSupportGeometry(&m_SectorSupportGeometry, &sectorSupport);
526  GearDir plane(gd);
527  plane.append("/Plane");
528  readPositionData(&m_PlanePosition, &plane);
529  readSizeData(&m_PlanePosition, &plane);
530  GearDir plasticSheet(gd);
531  plasticSheet.append("/PlasticSheet");
532  m_PlasticSheetGeometry.setWidth(plasticSheet.getLength("Width") * CLHEP::cm);
533  m_PlasticSheetGeometry.setDeltaL(plasticSheet.getLength("DeltaL") *
534  CLHEP::cm);
535  GearDir segmentSupport(gd);
536  segmentSupport.append("/SegmentSupport");
537  m_SegmentSupportGeometry.setTopWidth(
538  segmentSupport.getLength("TopWidth") * CLHEP::cm);
539  m_SegmentSupportGeometry.setTopThickness(
540  segmentSupport.getLength("TopThickness") * CLHEP::cm);
541  m_SegmentSupportGeometry.setMiddleWidth(
542  segmentSupport.getLength("MiddleWidth") * CLHEP::cm);
543  m_SegmentSupportGeometry.setMiddleThickness(
544  segmentSupport.getLength("MiddleThickness") * CLHEP::cm);
545  try {
546  m_SegmentSupportPosition =
547  new EKLMGeometry::SegmentSupportPosition[m_NSegmentSupportElementsSector];
548  } catch (std::bad_alloc& ba) {
549  B2FATAL(c_MemErr);
550  }
551  for (j = 0; j < m_NPlanes; j++) {
552  for (i = 0; i <= m_NSegments; i++) {
553  k = j * (m_NSegments + 1) + i;
554  GearDir segmentSupport2(segmentSupport);
555  name = "/SegmentSupportPlane[" + std::to_string(j + 1) + "]";
556  segmentSupport2.append(name.c_str());
557  name = "/SegmentSupport[" + std::to_string(i + 1) + "]";
558  segmentSupport2.append(name.c_str());
559  m_SegmentSupportPosition[k].setLength(
560  segmentSupport2.getLength("Length") * CLHEP::cm);
561  m_SegmentSupportPosition[k].setX(
562  segmentSupport2.getLength("X") * CLHEP::cm);
563  m_SegmentSupportPosition[k].setY(
564  segmentSupport2.getLength("Y") * CLHEP::cm);
565  m_SegmentSupportPosition[k].setZ(
566  segmentSupport2.getLength("Z") * CLHEP::cm);
567  m_SegmentSupportPosition[k].setDeltaLRight(
568  segmentSupport2.getLength("DeltaLRight") * CLHEP::cm);
569  m_SegmentSupportPosition[k].setDeltaLLeft(
570  segmentSupport2.getLength("DeltaLLeft") * CLHEP::cm);
571  }
572  }
573  readXMLDataStrips(gd);
574  GearDir shield(gd);
575  shield.append("/Shield");
576  m_ShieldGeometry.setThickness(shield.getLength("Thickness") * CLHEP::cm);
577  GearDir shieldDetailA(shield);
578  shieldDetailA.append("/Detail[@id=\"A\"]");
579  readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailA);
580  m_ShieldGeometry.setDetailA(shieldDetailGeometry);
581  GearDir shieldDetailB(shield);
582  shieldDetailB.append("/Detail[@id=\"B\"]");
583  readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailB);
584  m_ShieldGeometry.setDetailB(shieldDetailGeometry);
585  GearDir shieldDetailC(shield);
586  shieldDetailC.append("/Detail[@id=\"C\"]");
587  readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailC);
588  m_ShieldGeometry.setDetailC(shieldDetailGeometry);
589  GearDir shieldDetailD(shield);
590  shieldDetailD.append("/Detail[@id=\"D\"]");
591  readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailD);
592  m_ShieldGeometry.setDetailD(shieldDetailGeometry);
593  m_Geometry = new EKLMGeometry(*this);
594 }
595 
596 void EKLM::GeometryData::initializeFromDatabase()
597 {
598  DBObjPtr<EKLMGeometry> eklmGeometry;
599  if (!eklmGeometry.isValid())
600  B2FATAL("No EKLM geometry data in the database.");
601  EKLMGeometry::operator=(*eklmGeometry);
602  m_Geometry = new EKLMGeometry(*this);
603 }
604 
605 EKLM::GeometryData::GeometryData(enum DataSource dataSource,
606  const GearDir* gearDir)
607 {
608  m_Geometry = nullptr;
609  switch (dataSource) {
610  case c_Gearbox:
611  initializeFromGearbox(gearDir);
612  break;
613  case c_Database:
614  initializeFromDatabase();
615  break;
616  }
617  m_MinZForward = m_SolenoidZ + m_SectionPosition.getZ() -
618  0.5 * m_SectionPosition.getLength();
619  m_MaxZBackward = m_SolenoidZ - m_SectionPosition.getZ() +
620  0.5 * m_SectionPosition.getLength();
621  calculateSectorSupportGeometry();
622  fillStripIndexArrays();
623  calculateShieldGeometry();
624 }
625 
626 EKLM::GeometryData::~GeometryData()
627 {
628  if (m_Geometry != nullptr)
629  delete m_Geometry;
630  free(m_StripLenToAll);
631  free(m_StripAllToLen);
632 }
633 
634 void EKLM::GeometryData::saveToDatabase(const IntervalOfValidity& iov) const
635 {
636  Database::Instance().storeData("EKLMGeometry", m_Geometry, iov);
637 }
638 
639 bool EKLM::GeometryData::hitInEKLM(double z) const
640 {
641  double zMm;
642  zMm = z / Unit::cm * CLHEP::cm;
643  return (zMm > m_MinZForward) || (zMm < m_MaxZBackward);
644 }
645 
646 /*
647  * Note that numbers of elements are 0-based for all transformation functions.
648  */
649 void
650 EKLM::GeometryData::getSectionTransform(HepGeom::Transform3D* t, int n) const
651 {
652  if (n == 0)
653  *t = HepGeom::Translate3D(m_SectionPosition.getX(), m_SectionPosition.getY(),
654  -m_SectionPosition.getZ() + m_SolenoidZ);
655  else
656  *t = HepGeom::Translate3D(m_SectionPosition.getX(), m_SectionPosition.getY(),
657  m_SectionPosition.getZ() + m_SolenoidZ) *
658  HepGeom::RotateY3D(180.*CLHEP::deg);
659 }
660 
661 void
662 EKLM::GeometryData::getLayerTransform(HepGeom::Transform3D* t, int n) const
663 {
664  *t = HepGeom::Translate3D(0.0, 0.0, m_SectionPosition.getLength() / 2.0 -
665  (n + 1) * m_LayerShiftZ +
666  0.5 * m_LayerPosition.getLength());
667 }
668 
669 void
670 EKLM::GeometryData::getSectorTransform(HepGeom::Transform3D* t, int n) const
671 {
672  switch (n) {
673  case 0:
674  *t = HepGeom::Translate3D(0., 0., 0.);
675  break;
676  case 1:
677  *t = HepGeom::RotateY3D(180.0 * CLHEP::deg);
678  break;
679  case 2:
680  *t = HepGeom::RotateZ3D(90.0 * CLHEP::deg) *
681  HepGeom::RotateY3D(180.0 * CLHEP::deg);
682  break;
683  case 3:
684  *t = HepGeom::RotateZ3D(-90.0 * CLHEP::deg);
685  break;
686  }
687 }
688 
689 void
690 EKLM::GeometryData::getPlaneTransform(HepGeom::Transform3D* t, int n) const
691 {
692  if (n == 0)
693  *t = HepGeom::Translate3D(m_PlanePosition.getX(), m_PlanePosition.getY(),
694  m_PlanePosition.getZ()) *
695  HepGeom::Rotate3D(180. * CLHEP::deg,
696  HepGeom::Vector3D<double>(1., 1., 0.));
697  else
698  *t = HepGeom::Translate3D(m_PlanePosition.getX(), m_PlanePosition.getY(),
699  -m_PlanePosition.getZ());
700 }
701 
702 void
703 EKLM::GeometryData::getStripTransform(HepGeom::Transform3D* t, int n) const
704 {
705  *t = HepGeom::Translate3D(m_StripPosition[n].getX(),
706  m_StripPosition[n].getY(), 0.0);
707 }
708 
709 void
710 EKLM::GeometryData::getSheetTransform(HepGeom::Transform3D* t, int n) const
711 {
712  double y;
713  y = m_StripPosition[n].getY();
714  if (n % m_ElementNumbers->getNStripsSegment() == 0)
715  y = y + 0.5 * m_PlasticSheetGeometry.getDeltaL();
716  else if (n % m_ElementNumbers->getNStripsSegment() ==
717  m_ElementNumbers->getNStripsSegment() - 1)
718  y = y - 0.5 * m_PlasticSheetGeometry.getDeltaL();
719  *t = HepGeom::Translate3D(m_StripPosition[n].getX(), y, 0.0);
720 }
721 
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:35
Belle2::EKLM::GeometryData::~GeometryData
~GeometryData()
Destructor.
Definition: GeometryData.cc:626
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner4Z
void setCorner4Z(double corner4Z)
Set corner 4 Z coordinate.
Definition: EKLMGeometry.h:576
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::EKLMGeometry
Class to store EKLM geometry data in the database.
Definition: EKLMGeometry.h:39
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner3Z
void setCorner3Z(double corner3Z)
Set corner 3 Z coordinate.
Definition: EKLMGeometry.h:508
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner4LX
void setCorner4LX(double corner4LX)
Set corner 4 X length.
Definition: EKLMGeometry.h:525
Belle2::EKLMGeometry::ShieldDetailGeometry::setPoint
void setPoint(int i, const Point &point)
Set point.
Definition: EKLMGeometry.cc:142
Belle2::EKLM::GeometryData::initializeFromDatabase
void initializeFromDatabase()
Initialize from database.
Definition: GeometryData.cc:596
Belle2::gearbox::Interface::getInt
int getInt(const std::string &path="") const noexcept(false)
Get the parameter path as a int.
Definition: Interface.cc:70
Belle2::EKLM::GeometryData::readEndcapStructureGeometry
void readEndcapStructureGeometry(const GearDir &gd)
Read section structure geometry data.
Definition: GeometryData.cc:471
Belle2::EKLM::GeometryData::fillStripIndexArrays
void fillStripIndexArrays()
Fill strip index arrays.
Definition: GeometryData.cc:204
Belle2::EKLM::GeometryData::getLayerTransform
void getLayerTransform(HepGeom::Transform3D *t, int n) const
Get layer transformation.
Definition: GeometryData.cc:662
Belle2::Database::storeData
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition: Database.cc:152
Belle2::EKLMGeometry::SectorSupportGeometry
Sector support geometry data.
Definition: EKLMGeometry.h:249
Belle2::EKLMGeometry::Point
2D point.
Definition: EKLMGeometry.h:1326
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner2Z
void setCorner2Z(double corner2Z)
Set corner 2 Z coordinate.
Definition: EKLMGeometry.h:440
Belle2::EKLMGeometry::ShieldDetailGeometry
Shield layer detail geometry data.
Definition: EKLMGeometry.h:1385
Belle2::EKLMGeometry::ElementPosition
Position information for the elements of detector.
Definition: EKLMGeometry.h:110
Belle2::EKLM::GeometryData::calculateSectorSupportGeometry
void calculateSectorSupportGeometry()
Calculate sector support geometry data.
Definition: GeometryData.cc:119
Belle2::EKLMGeometry::ShieldDetailGeometry::setLengthX
void setLengthX(double lengthX)
Set X length.
Definition: EKLMGeometry.h:1421
Belle2::EKLM::Line2D
2D line.
Definition: Line2D.h:40
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner3LX
void setCorner3LX(double corner3LX)
Set corner 3 X length.
Definition: EKLMGeometry.h:457
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner1Z
void setCorner1Z(double corner1Z)
Set corner 1 Z coordinate.
Definition: EKLMGeometry.h:372
Belle2::EKLM::Circle2D
2D circle.
Definition: Circle2D.h:35
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner1Thickness
void setCorner1Thickness(double corner1Thickness)
Set corner 1 thickness.
Definition: EKLMGeometry.h:355
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::EKLMGeometry::operator=
EKLMGeometry & operator=(const EKLMGeometry &geometry)
Operator =.
Definition: EKLMGeometry.cc:244
Belle2::EKLMGeometry::SegmentSupportPosition
Segment support position.
Definition: EKLMGeometry.h:1048
Belle2::EKLM::GeometryData::getSheetTransform
void getSheetTransform(HepGeom::Transform3D *t, int n) const
Get plastic sheet element transformation.
Definition: GeometryData.cc:710
Belle2::EKLMGeometry::SectorSupportGeometry::setThickness
void setThickness(double thickness)
Set thickness.
Definition: EKLMGeometry.h:270
Belle2::EKLMGeometry::ShieldDetailGeometry::getLengthY
double getLengthY() const
Get Y length.
Definition: EKLMGeometry.h:1429
Belle2::EKLM::GeometryData::getPlaneTransform
void getPlaneTransform(HepGeom::Transform3D *t, int n) const
Get plane transformation.
Definition: GeometryData.cc:690
Belle2::EKLMGeometry::SectorSupportGeometry::setDeltaLY
void setDeltaLY(double deltaLY)
Set outerR - Y of upper edge of BoxY.
Definition: EKLMGeometry.h:287
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner3LY
void setCorner3LY(double corner3LY)
Set corner 3 Y length.
Definition: EKLMGeometry.h:474
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::EKLM::GeometryData::saveToDatabase
void saveToDatabase(const IntervalOfValidity &iov) const
Save geometry data to database.
Definition: GeometryData.cc:634
Belle2::EKLMGeometry::ShieldDetailGeometry::getPoint
const Point * getPoint(int i) const
Get point.
Definition: EKLMGeometry.cc:135
Belle2::EKLMGeometry::ElementPosition::setOuterR
void setOuterR(double outerR)
Set outer radius.
Definition: EKLMGeometry.h:148
Belle2::EKLM::GeometryData::hitInEKLM
bool hitInEKLM(double z) const
Check if z coordinate may be in EKLM.
Definition: GeometryData.cc:639
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner4Thickness
void setCorner4Thickness(double corner4Thickness)
Set corner 4 thickness.
Definition: EKLMGeometry.h:559
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::GearDir::append
void append(const std::string &path)
Append something to the current path, modifying the GearDir in place.
Definition: GearDir.h:62
Belle2::EKLMGeometry::Point::getX
double getX() const
Get X coordinate.
Definition: EKLMGeometry.h:1338
Belle2::EKLM::GeometryData
EKLM geometry data.
Definition: GeometryData.h:40
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner2Thickness
void setCorner2Thickness(double corner2Thickness)
Set corner 2 thickness.
Definition: EKLMGeometry.h:423
Belle2::EKLM::GeometryData::DataSource
DataSource
Geometry data source.
Definition: GeometryData.h:45
Belle2::EKLMGeometry::ShieldDetailGeometry::setLengthY
void setLengthY(double lengthY)
Set Y length.
Definition: EKLMGeometry.h:1438
prepareAsicCrosstalkSimDB.u
u
merged u1 and u2
Definition: prepareAsicCrosstalkSimDB.py:46
Belle2::EKLM::GeometryData::initializeFromGearbox
void initializeFromGearbox(const GearDir *gearDir)
Initialize from Gearbox (XML).
Definition: GeometryData.cc:479
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner1LX
void setCorner1LX(double corner1LX)
Set corner 1 X length.
Definition: EKLMGeometry.h:321
Belle2::EKLM::GeometryData::GeometryData
GeometryData(enum DataSource dataSource, const GearDir *gearDir)
Constructor.
Definition: GeometryData.cc:605
Belle2::gearbox::Interface::getLength
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:261
Belle2::EKLM::GeometryData::calculateShieldGeometry
void calculateShieldGeometry()
Calculate shield geometry data.
Definition: GeometryData.cc:387
Belle2::EKLMGeometry::ShieldDetailGeometry::setNPoints
void setNPoints(int nPoints)
Set number of points.
Definition: EKLMGeometry.cc:121
Belle2::DBAccessorBase::isValid
bool isValid() const
Check whether a valid object was obtained from the database.
Definition: DBAccessorBase.h:75
Belle2::EKLM::Line2D::findIntersection
int findIntersection(const Line2D &line, HepGeom::Point3D< double > *intersection) const
Find intersection with a line.
Definition: Line2D.cc:30
Belle2::EKLMGeometry::ShieldDetailGeometry::getLengthX
double getLengthX() const
Get X length.
Definition: EKLMGeometry.h:1412
Belle2::EKLM::GeometryData::getStripTransform
void getStripTransform(HepGeom::Transform3D *t, int n) const
Get strip transformation.
Definition: GeometryData.cc:703
Belle2::Database::Instance
static Database & Instance()
Instance of a singleton Database.
Definition: Database.cc:54
Belle2::EKLMGeometry::ElementPosition::setZ
void setZ(double z)
Set Z coordinate.
Definition: EKLMGeometry.h:216
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner2LY
void setCorner2LY(double corner2LY)
Set corner 2 Y length.
Definition: EKLMGeometry.h:406
Belle2::EKLMGeometry::ElementPosition::setX
void setX(double x)
Set X coordinate.
Definition: EKLMGeometry.h:182
Belle2::EKLMGeometry::Point::getY
double getY() const
Get Y coordinate.
Definition: EKLMGeometry.h:1355
HepGeom::Point3D< double >
Belle2::EKLMGeometry::ElementPosition::setLength
void setLength(double length)
Set length.
Definition: EKLMGeometry.h:165
Belle2::EKLMGeometry::ElementPosition::setInnerR
void setInnerR(double innerR)
Set inner radius.
Definition: EKLMGeometry.h:131
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner1Width
void setCorner1Width(double corner1Width)
Set corner 1 width.
Definition: EKLMGeometry.h:338
Belle2::EKLMGeometry::ElementPosition::setY
void setY(double y)
Set Y coordinate.
Definition: EKLMGeometry.h:199
Belle2::EKLM::GeometryData::getSectorTransform
void getSectorTransform(HepGeom::Transform3D *t, int n) const
Get sector transformation.
Definition: GeometryData.cc:670
Belle2::EKLM::GeometryData::readXMLDataStrips
void readXMLDataStrips(const GearDir &gd)
Read strip parameters from XML database.
Definition: GeometryData.cc:260
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner2LX
void setCorner2LX(double corner2LX)
Set corner 2 X length.
Definition: EKLMGeometry.h:389
Belle2::EKLMGeometry::SectorSupportGeometry::setCornerX
void setCornerX(double cornerX)
Set coordinate X of corner 1.
Definition: EKLMGeometry.h:304
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner4LY
void setCorner4LY(double corner4LY)
Set corner 4 Y length.
Definition: EKLMGeometry.h:542
Belle2::EKLM::GeometryData::getSectionTransform
void getSectionTransform(HepGeom::Transform3D *t, int n) const
Get section transformation.
Definition: GeometryData.cc:650
Belle2::EKLM::GeometryData::Instance
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
Definition: GeometryData.cc:35
Belle2::GearDir::getNumberNodes
virtual int getNumberNodes(const std::string &path="") const override
Return the number of nodes a given path will expand to.
Definition: GearDir.h:68
Belle2::EKLMGeometry::SectorSupportGeometry::setCorner3Thickness
void setCorner3Thickness(double corner3Thickness)
Set corner 3 thickness.
Definition: EKLMGeometry.h:491