10#include <klm/eklm/geometry/GeometryData.h>
13#include <klm/eklm/geometry/Circle2D.h>
14#include <klm/eklm/geometry/Line2D.h>
17#include <framework/database/DBObjPtr.h>
18#include <framework/database/Database.h>
19#include <framework/logging/Logger.h>
22#include <CLHEP/Geometry/Point3D.h>
23#include <CLHEP/Units/SystemOfUnits.h>
30static const char c_MemErr[] =
"Memory allocation error.";
68static void readSectorSupportGeometry(
97static void readShieldDetailGeometry(
98 EKLMGeometry::ShieldDetailGeometry* sdg, GearDir* gd)
102 EKLMGeometry::Point p;
107 for (i = 0; i < n; i++) {
109 name =
"/Point[" + std::to_string(i + 1) +
"]";
111 p.setX(point.getLength(
"X") * CLHEP::cm);
112 p.setY(point.getLength(
"Y") * CLHEP::cm);
179 p.setX(intersections[1].x());
192 p.setY(intersections[1].y());
197static bool compareLength(
double a,
double b)
204 const char err[] =
"Strip sorting algorithm error.";
207 std::vector<double> strips;
208 std::vector<double>::iterator it;
209 std::map<double, int> mapLengthStrip;
210 std::map<double, int> mapLengthStrip2;
211 std::map<double, int>::iterator itm;
214 mapLengthStrip.insert(
217 sort(strips.begin(), strips.end(), compareLength);
220 for (it = strips.begin(); it != strips.end(); ++it) {
231 itm = mapLengthStrip.find(l);
232 if (itm == mapLengthStrip.end())
235 mapLengthStrip2.insert(std::pair<double, int>(l, i));
236 for (it = strips.begin(); it != strips.end(); ++it) {
240 itm = mapLengthStrip.find(l);
241 if (itm == mapLengthStrip.end())
244 mapLengthStrip2.insert(std::pair<double, int>(l, i));
252 if (itm == mapLengthStrip2.end())
269 Strips.
getLength(
"NoScintillationThickness") * CLHEP::cm);
273 }
catch (std::bad_alloc& ba) {
278 name =
"/Strip[" + std::to_string(i + 1) +
"]";
279 StripContent.
append(name);
298 double r,
double kx,
double ky,
299 double& dx,
double& dy)
304 double a, b, c, d, t, maxt = 0, x1, y1, x2, y2, u;
311 a = kx * kx + ky * ky;
312 intersection =
false;
313 for (i = 0; i < nPoints; i++) {
316 b = 2.0 * (kx * x1 + ky * y1);
317 c = x1 * x1 + y1 * y1 - r * r;
318 d = b * b - 4.0 * a * c;
320 t = (-b +
sqrt(d)) / (2.0 * a);
331 B2FATAL(
"Shield layer geometry calculation error.");
342 for (i = 0; i < nPoints; i++) {
345 if (i < nPoints - 1) {
346 x2 = points[i + 1].x();
347 y2 = points[i + 1].y();
352 a = (x2 - x1) * ky - (y2 - y1) * kx;
355 b = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
356 t = (x1 * y2 - x2 * y1 + r *
sqrt(b)) / a;
362 u = -((x2 - x1) * (x1 + kx * t) + (y2 - y1) * (y2 + ky * t)) / b;
377static void EKLMPointToCLHEP(
const EKLMGeometry::Point* pointEKLM,
378 HepGeom::Point3D<double>& pointCLHEP)
380 pointCLHEP.setX(pointEKLM->
getX());
381 pointCLHEP.setY(pointEKLM->
getY());
388 double r, l, dx, dy, xCenter, yCenter;
389 const double asqrt2 = 1.0 /
sqrt(2.0);
398 EKLMPointToCLHEP(detailA->
getPoint(0), points[0]);
399 EKLMPointToCLHEP(detailA->
getPoint(1), points[1]);
400 EKLMPointToCLHEP(detailA->
getPoint(2), points[2]);
401 EKLMPointToCLHEP(detailA->
getPoint(3), points[3]);
402 EKLMPointToCLHEP(detailA->
getPoint(4), points[4]);
406 EKLMPointToCLHEP(detailA->
getPoint(5), points[6]);
407 EKLMPointToCLHEP(detailA->
getPoint(6), points[7]);
409 xCenter = -asqrt2 * l;
410 yCenter = asqrt2 * l;
411 for (i = 0; i < 8; i++)
412 points[i] = HepGeom::Translate3D(xCenter, yCenter, 0) *
413 HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
414 HepGeom::Translate3D(-detailA->
getLengthX() / 2,
417 getDetailDxDy(points, 8, r, 1, 1, dx, dy);
432 EKLMPointToCLHEP(detailB->
getPoint(0), points[4]);
433 EKLMPointToCLHEP(detailB->
getPoint(1), points[5]);
434 EKLMPointToCLHEP(detailB->
getPoint(2), points[6]);
435 EKLMPointToCLHEP(detailB->
getPoint(3), points[7]);
437 for (i = 0; i < 8; i++)
438 points[i] = HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
439 HepGeom::Translate3D(-detailB->
getLengthX() / 2,
442 getDetailDxDy(points, 8, r, 1, 1, dx, dy);
445 EKLMPointToCLHEP(detailC->
getPoint(0), points[0]);
446 EKLMPointToCLHEP(detailC->
getPoint(1), points[1]);
447 EKLMPointToCLHEP(detailC->
getPoint(2), points[2]);
448 EKLMPointToCLHEP(detailC->
getPoint(3), points[3]);
449 EKLMPointToCLHEP(detailC->
getPoint(4), points[4]);
453 EKLMPointToCLHEP(detailC->
getPoint(5), points[6]);
454 EKLMPointToCLHEP(detailC->
getPoint(6), points[7]);
456 xCenter = asqrt2 * l;
457 yCenter = -asqrt2 * l;
458 for (i = 0; i < 8; i++)
459 points[i] = HepGeom::Translate3D(xCenter, yCenter, 0) *
460 HepGeom::RotateZ3D(-45.0 * CLHEP::deg) *
461 HepGeom::RotateY3D(180.0 * CLHEP::deg) *
462 HepGeom::Translate3D(-detailC->
getLengthX() / 2,
465 getDetailDxDy(points, 8, r, 1, 1, dx, dy);
472 d.append(
"/EndcapStructure");
511 section.append(
"/Section");
515 layer.append(
"/Layer");
519 sector.append(
"/Sector");
522 sectorSupport.
append(
"/SectorSupport");
527 plane.append(
"/Plane");
531 plasticSheet.
append(
"/PlasticSheet");
536 segmentSupport.
append(
"/SegmentSupport");
538 segmentSupport.
getLength(
"TopWidth") * CLHEP::cm);
540 segmentSupport.
getLength(
"TopThickness") * CLHEP::cm);
542 segmentSupport.
getLength(
"MiddleWidth") * CLHEP::cm);
544 segmentSupport.
getLength(
"MiddleThickness") * CLHEP::cm);
548 }
catch (std::bad_alloc& ba) {
554 GearDir segmentSupport2(segmentSupport);
555 name =
"/SegmentSupportPlane[" + std::to_string(j + 1) +
"]";
556 segmentSupport2.append(name);
557 name =
"/SegmentSupport[" + std::to_string(i + 1) +
"]";
558 segmentSupport2.append(name);
560 segmentSupport2.getLength(
"Length") * CLHEP::cm);
562 segmentSupport2.getLength(
"X") * CLHEP::cm);
564 segmentSupport2.getLength(
"Y") * CLHEP::cm);
566 segmentSupport2.getLength(
"Z") * CLHEP::cm);
568 segmentSupport2.getLength(
"DeltaLRight") * CLHEP::cm);
570 segmentSupport2.getLength(
"DeltaLLeft") * CLHEP::cm);
578 shieldDetailA.
append(
"/Detail[@id=\"A\"]");
579 readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailA);
582 shieldDetailB.
append(
"/Detail[@id=\"B\"]");
583 readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailB);
586 shieldDetailC.
append(
"/Detail[@id=\"C\"]");
587 readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailC);
590 shieldDetailD.
append(
"/Detail[@id=\"D\"]");
591 readShieldDetailGeometry(&shieldDetailGeometry, &shieldDetailD);
600 B2FATAL(
"No EKLM geometry data in the database.");
609 switch (dataSource) {
658 HepGeom::RotateY3D(180.*CLHEP::deg);
674 *t = HepGeom::Translate3D(0., 0., 0.);
677 *t = HepGeom::RotateY3D(180.0 * CLHEP::deg);
680 *t = HepGeom::RotateZ3D(90.0 * CLHEP::deg) *
681 HepGeom::RotateY3D(180.0 * CLHEP::deg);
684 *t = HepGeom::RotateZ3D(-90.0 * CLHEP::deg);
695 HepGeom::Rotate3D(180. * CLHEP::deg,
696 HepGeom::Vector3D<double>(1., 1., 0.));
bool isValid() const
Check whether a valid object was obtained from the database.
Class for accessing objects in the database.
Position information for the elements of detector.
void setInnerR(double innerR)
Set inner radius.
void setOuterR(double outerR)
Set outer radius.
void setZ(double z)
Set Z coordinate.
void setLength(double length)
Set length.
void setY(double y)
Set Y coordinate.
void setX(double x)
Set X coordinate.
double getX() const
Get X coordinate.
double getY() const
Get Y coordinate.
Sector support geometry data.
void setCornerX(double cornerX)
Set coordinate X of corner 1.
void setCorner1LX(double corner1LX)
Set corner 1 X length.
void setCorner2LY(double corner2LY)
Set corner 2 Y length.
void setCorner4Thickness(double corner4Thickness)
Set corner 4 thickness.
void setThickness(double thickness)
Set thickness.
void setCorner1Z(double corner1Z)
Set corner 1 Z coordinate.
void setCorner1Thickness(double corner1Thickness)
Set corner 1 thickness.
void setCorner2Thickness(double corner2Thickness)
Set corner 2 thickness.
void setCorner3LY(double corner3LY)
Set corner 3 Y length.
void setDeltaLY(double deltaLY)
Set outerR - Y of upper edge of BoxY.
void setCorner4LY(double corner4LY)
Set corner 4 Y length.
void setCorner3LX(double corner3LX)
Set corner 3 X length.
void setCorner4Z(double corner4Z)
Set corner 4 Z coordinate.
void setCorner3Thickness(double corner3Thickness)
Set corner 3 thickness.
void setCorner3Z(double corner3Z)
Set corner 3 Z coordinate.
void setCorner2LX(double corner2LX)
Set corner 2 X length.
void setCorner2Z(double corner2Z)
Set corner 2 Z coordinate.
void setCorner1Width(double corner1Width)
Set corner 1 width.
void setCorner4LX(double corner4LX)
Set corner 4 X length.
Segment support position.
Shield layer detail geometry data.
void setPoint(int i, const Point &point)
Set point.
void setLengthY(double lengthY)
Set Y length.
const Point * getPoint(int i) const
Get point.
double getLengthY() const
Get Y length.
double getLengthX() const
Get X length.
void setLengthX(double lengthX)
Set X length.
void setNPoints(int nPoints)
Set number of points.
int m_NSegments
Number of segments in one plane.
EKLMGeometry()
Constructor.
ShieldGeometry m_ShieldGeometry
Shield layer details geometry data.
int m_NLayers
Number of layers in one section.
bool m_BeamBackgroundStudy
ROOT streamer.
SectorSupportGeometry m_SectorSupportGeometry
Sector support geometry data.
ElementPosition m_SectionPosition
Position data for sections.
int m_NPlanes
Number of planes in one sector.
ElementPosition m_SectorSupportPosition
Position data for sector support structure.
double m_SolenoidZ
Solenoid center Z coordinate.
ElementPosition m_PlanePosition
Position data for planes.
ElementPosition m_LayerPosition
Position data for layers.
int m_NStrips
Number of strips in one plane.
double m_LayerShiftZ
Z distance between two layers.
int m_NSectors
Number of sectors in one layer.
const EKLMElementNumbers * m_ElementNumbers
Element numbers.
int * m_NDetectorLayers
Number of detector layers.
int m_NSections
Number of sections.
PlasticSheetGeometry m_PlasticSheetGeometry
Plastic sheet geometry data.
ElementPosition m_SectorPosition
Position data for sectors.
SegmentSupportPosition * m_SegmentSupportPosition
Position data for segment support structure.
ElementPosition * m_StripPosition
Position data for strips.
StripGeometry m_StripGeometry
Strip geometry data.
EndcapStructureGeometry m_EndcapStructureGeometry
Section structure geometry data.
EKLMGeometry & operator=(const EKLMGeometry &geometry)
Operator =.
SegmentSupportGeometry m_SegmentSupportGeometry
Segment support geometry data.
void checkDetectorLayerNumber(int section, int layer) const
Check if number of detector layers is correct (fatal error if not).
int m_NSegmentSupportElementsSector
Number of segment support elements in one sector.
void getSectorTransform(HepGeom::Transform3D *t, int n) const
Get sector transformation.
int * m_StripLenToAll
Number of strip in position-based array.
void fillStripIndexArrays()
Fill strip index arrays.
void getSheetTransform(HepGeom::Transform3D *t, int n) const
Get plastic sheet element transformation.
bool hitInEKLM(double z) const
Check if z coordinate may be in EKLM.
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
void calculateSectorSupportGeometry()
Calculate sector support geometry data.
void saveToDatabase(const IntervalOfValidity &iov) const
Save geometry data to database.
GeometryData(enum DataSource dataSource, const GearDir *gearDir)
Constructor.
void initializeFromGearbox(const GearDir *gearDir)
Initialize from Gearbox (XML).
void getStripTransform(HepGeom::Transform3D *t, int n) const
Get strip transformation.
void readEndcapStructureGeometry(const GearDir &gd)
Read section structure geometry data.
void getSectionTransform(HepGeom::Transform3D *t, int n) const
Get section transformation.
void getLayerTransform(HepGeom::Transform3D *t, int n) const
Get layer transformation.
void calculateShieldGeometry()
Calculate shield geometry data.
double m_MaxZBackward
Maximal z coordinate of the backward section.
EKLMGeometry * m_Geometry
Copy of data in this class used to write it to database.
double m_MinZForward
Minimal z coordinate of the forward section.
void initializeFromDatabase()
Initialize from database.
DataSource
Geometry data source.
@ c_Gearbox
Gearbox (XML).
int m_nStripDifferent
Number of strips with different lengths in one plane.
void getPlaneTransform(HepGeom::Transform3D *t, int n) const
Get plane transformation.
~GeometryData()
Destructor.
int * m_StripAllToLen
Number of strip in length-based array.
void readXMLDataStrips(const GearDir &gd)
Read strip parameters from XML database.
int findIntersection(const Line2D &line, HepGeom::Point3D< double > *intersection) const
Find intersection with a line.
GearDir is the basic class used for accessing the parameter store.
void append(const std::string &path)
Append something to the current path, modifying the GearDir in place.
virtual int getNumberNodes(const std::string &path="") const override
Return the number of nodes a given path will expand to.
A class that describes the interval of experiments/runs for which an object in the database is valid.
static const double cm
Standard units with the value = 1.
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
bool getBool(const std::string &path="") const noexcept(false)
Get the parameter path as a bool.
int getInt(const std::string &path="") const noexcept(false)
Get the parameter path as a int.
static Database & Instance()
Instance of a singleton Database.
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
double tan(double a)
tan for double
double sqrt(double a)
sqrt for double
Abstract base class for different kinds of events.