28 int iSection, iLayer, iSector, iPlane, iSegment, iStrip, sector, segment;
29 int nSections, nLayers, nSectors, nPlanes, nStrips, nSegments, nStripsSegment;
34 nSections =
m_GeoDat->getNSections();
39 nSegments =
m_GeoDat->getNSegments();
41 m_Section =
new HepGeom::Transform3D[nSections];
42 m_Layer =
new HepGeom::Transform3D*[nSections];
43 m_Sector =
new HepGeom::Transform3D** [nSections];
44 m_Plane =
new HepGeom::Transform3D** *[nSections];
46 m_Segment =
new HepGeom::Transform3D**** [nSections];
47 m_Strip =
new HepGeom::Transform3D**** [nSections];
49 for (iSection = 0; iSection < nSections; iSection++) {
51 nDetectorLayers =
m_GeoDat->getNDetectorLayers(iSection + 1);
52 m_Layer[iSection] =
new HepGeom::Transform3D[nLayers];
53 m_Sector[iSection] =
new HepGeom::Transform3D*[nLayers];
54 m_Plane[iSection] =
new HepGeom::Transform3D** [nLayers];
56 m_Segment[iSection] =
new HepGeom::Transform3D** *[nLayers];
57 m_Strip[iSection] =
new HepGeom::Transform3D** *[nLayers];
59 for (iLayer = 0; iLayer < nLayers; iLayer++) {
61 m_Sector[iSection][iLayer] =
new HepGeom::Transform3D[nSectors];
62 if (iLayer < nDetectorLayers) {
63 m_Plane[iSection][iLayer] =
new HepGeom::Transform3D*[nSectors];
65 new HepGeom::Transform3D*[nSectors];
66 m_Segment[iSection][iLayer] =
new HepGeom::Transform3D** [nSectors];
67 m_Strip[iSection][iLayer] =
new HepGeom::Transform3D** [nSectors];
68 m_StripInverse[iSection][iLayer] =
new HepGeom::Transform3D** [nSectors];
70 for (iSector = 0; iSector < nSectors; iSector++) {
73 if (iLayer >= nDetectorLayers)
75 m_Plane[iSection][iLayer][iSector] =
new HepGeom::Transform3D[nPlanes];
77 new HepGeom::Transform3D[nPlanes];
79 new HepGeom::Transform3D*[nPlanes];
80 m_Strip[iSection][iLayer][iSector] =
new HepGeom::Transform3D*[nPlanes];
82 new HepGeom::Transform3D*[nPlanes];
83 for (iPlane = 0; iPlane < nPlanes; iPlane++) {
85 &
m_Plane[iSection][iLayer][iSector][iPlane], iPlane);
87 HepGeom::Translate3D(0, 0, 0);
88 m_Segment[iSection][iLayer][iSector][iPlane] =
89 new HepGeom::Transform3D[nSegments];
90 for (iSegment = 0; iSegment < nSegments; iSegment++) {
91 m_Segment[iSection][iLayer][iSector][iPlane][iSegment] =
92 HepGeom::Translate3D(0, 0, 0);
94 m_Strip[iSection][iLayer][iSector][iPlane] =
95 new HepGeom::Transform3D[nStrips];
97 new HepGeom::Transform3D[nStrips];
98 for (iStrip = 0; iStrip < nStrips; iStrip++) {
100 &
m_Strip[iSection][iLayer][iSector][iPlane][iStrip], iStrip);
107 if (displacementType !=
c_None) {
108 std::string payload, segmentPayload;
110 payload =
"EKLMDisplacement";
111 segmentPayload =
"EKLMSegmentDisplacement";
113 payload =
"EKLMAlignment";
114 segmentPayload =
"EKLMSegmentAlignment";
118 if (!alignment.isValid())
119 B2FATAL(
"No EKLM displacement (alignment) data.");
121 if (!alignmentChecker.
checkAlignment(&(*alignment), &(*segmentAlignment)))
122 B2FATAL(
"EKLM displacement data are incorrect, overlaps exist.");
124 for (iSection = 1; iSection <= nSections; iSection++) {
125 nDetectorLayers =
m_GeoDat->getNDetectorLayers(iSection);
126 for (iLayer = 1; iLayer <= nDetectorLayers; iLayer++) {
127 for (iSector = 1; iSector <= nSectors; iSector++) {
130 alignment->getModuleAlignment(sector);
131 if (sectorAlignment ==
nullptr)
132 B2FATAL(
"Incomplete EKLM displacement (alignment) data.");
133 for (iPlane = 1; iPlane <= nPlanes; iPlane++) {
140 HepGeom::Translate3D(
147 B2FATAL(
"Missing m_PlaneDisplacement allocation for section/layer/sector in TransformData.cc");
151 HepGeom::Translate3D(
157 for (iSegment = 1; iSegment <= nSegments; iSegment++) {
159 iSection, iLayer, iSector, iPlane, iSegment);
161 segmentAlignment->getSegmentAlignment(segment);
162 if (segmentAlignmentData ==
nullptr)
163 B2FATAL(
"Incomplete EKLM displacement (alignment) data.");
164 m_Segment[iSection - 1][iLayer - 1][iSector - 1][iPlane - 1]
166 HepGeom::Translate3D(
169 m_Segment[iSection - 1][iLayer - 1][iSector - 1][iPlane - 1]
173 for (iStrip = 1; iStrip <= nStripsSegment; iStrip++) {
174 m_Strip[iSection - 1][iLayer - 1][iSector - 1][iPlane - 1]
175 [nStripsSegment * (iSegment - 1) + iStrip - 1] =
176 HepGeom::Translate3D(
180 m_Strip[iSection - 1][iLayer - 1][iSector - 1][iPlane - 1]
181 [nStripsSegment * (iSegment - 1) + iStrip - 1] *
197 int iSection, iLayer, iSector, iPlane;
199 int nSections, nLayers, nDetectorLayers, nSectors, nPlanes;
200 nSections =
m_GeoDat->getNSections();
204 for (iSection = 0; iSection < nSections; iSection++) {
205 nDetectorLayers =
m_GeoDat->getNDetectorLayers(iSection + 1);
206 for (iLayer = 0; iLayer < nLayers; iLayer++) {
207 delete[]
m_Sector[iSection][iLayer];
208 if (iLayer >= nDetectorLayers)
210 for (iSector = 0; iSector < nSectors; iSector++) {
211 for (iPlane = 0; iPlane < nPlanes; iPlane++) {
212 delete[]
m_Segment[iSection][iLayer][iSector][iPlane];
213 delete[]
m_Strip[iSection][iLayer][iSector][iPlane];
216 delete[]
m_Plane[iSection][iLayer][iSector];
218 delete[]
m_Segment[iSection][iLayer][iSector];
219 delete[]
m_Strip[iSection][iLayer][iSector];
222 delete[]
m_Plane[iSection][iLayer];
225 delete[]
m_Strip[iSection][iLayer];
248 int iSection, iLayer, iSector, iPlane, iSegment, iStrip;
250 int nSections, nLayers, nDetectorLayers, nSectors, nPlanes, nSegments, nStrips;
251 nSections =
m_GeoDat->getNSections();
255 nSegments =
m_GeoDat->getNSegments();
257 for (iSection = 0; iSection < nSections; iSection++) {
258 nDetectorLayers =
m_GeoDat->getNDetectorLayers(iSection + 1);
259 for (iLayer = 0; iLayer < nLayers; iLayer++) {
261 for (iSector = 0; iSector < nSectors; iSector++) {
262 m_Sector[iSection][iLayer][iSector] =
264 if (iLayer >= nDetectorLayers)
266 for (iPlane = 0; iPlane < nPlanes; iPlane++) {
269 m_Plane[iSection][iLayer][iSector]) {
270 m_Plane[iSection][iLayer][iSector][iPlane] =
271 m_Sector[iSection][iLayer][iSector] *
272 m_Plane[iSection][iLayer][iSector][iPlane];
275 B2FATAL(
"Missing m_Plane allocation for section/layer/sector");
277 for (iSegment = 0; iSegment < nSegments; iSegment++) {
278 m_Segment[iSection][iLayer][iSector][iPlane][iSegment] =
279 m_Plane[iSection][iLayer][iSector][iPlane] *
281 m_Segment[iSection][iLayer][iSector][iPlane][iSegment];
283 for (iStrip = 0; iStrip < nStrips; iStrip++) {
284 m_Strip[iSection][iLayer][iSector][iPlane][iStrip] =
285 m_Plane[iSection][iLayer][iSector][iPlane] *
287 m_Strip[iSection][iLayer][iSector][iPlane][iStrip];
289 m_Strip[iSection][iLayer][iSector][iPlane][iStrip].inverse();
364 double* d1,
double* d2,
double* sd,
400 HepGeom::Vector3D<double> v1 = s1_2g - s1_1g;
401 HepGeom::Vector3D<double> v2 = s2_2g - s2_1g;
402 HepGeom::Vector3D<double> d = s1_1g - s2_1g;
403 double v1sq = v1.mag2();
404 double v2sq = v2.mag2();
405 double v1dv2 = v1.dot(v2);
406 double ddv1 = d.dot(v1);
407 double ddv2 = d.dot(v2);
408 double den = v1sq * v2sq - v1dv2 * v1dv2;
409 double t1 = (v1dv2 * ddv2 - v2sq * ddv1) / den;
410 double t2 = (- v1dv2 * ddv1 + v1sq * ddv2) / den;
413 if (t1 < 0.0 || t1 > 1.0)
415 if (t2 < 0.0 || t2 > 1.0)
421 *d1 = s1_2g.distance(s1_cg) / CLHEP::mm *
Unit::mm;
422 *d2 = s2_2g.distance(s2_cg) / CLHEP::mm *
Unit::mm;
423 *cross = 0.5 * (s1_cg + s2_cg) / CLHEP::mm *
Unit::mm;
424 *sd = s1_cg.distance(s2_cg) / CLHEP::mm *
Unit::mm;
425 if (s2_cg.mag2() < s1_cg.mag2())
458 int section, layer, sector, plane, segment, strip, stripSegment, stripGlobal;
460 int nLayers, nPlanes, nSegments, nStripsSegment, minDistanceSegment;
461 double solenoidCenter, firstLayerCenter, layerShift;
463 double x, y, z, l, minY, maxY;
464 double minDistance = 0, minDistanceNew, stripWidth;
473 (
m_GeoDat->getSectionPosition()->getZ()
474 - 0.5 *
m_GeoDat->getSectionPosition()->getLength()
476 - 0.5 *
m_GeoDat->getLayerPosition()->getLength()) /
480 layer = round((z - firstLayerCenter) / layerShift) + 1;
483 nLayers =
m_GeoDat->getNDetectorLayers(section);
488 nSegments =
m_GeoDat->getNSegments();
490 stripWidth =
m_GeoDat->getStripGeometry()->getWidth() / CLHEP::cm *
Unit::cm;
491 minY = -stripWidth / 2;
492 maxY = (double(nStripsSegment) - 0.5) * stripWidth;
493 for (plane = 1; plane <= nPlanes; plane++) {
494 minDistanceSegment = 1;
495 for (segment = 1; segment <= nSegments; segment++) {
496 strip = (segment - 1) * nStripsSegment;
498 [sector - 1][plane - 1][strip] * intersectionClhep;
499 y = intersectionLocal.y() / CLHEP::cm *
Unit::cm;
501 minDistanceNew = minY - y;
502 }
else if (y > maxY) {
503 minDistanceNew = y - maxY;
506 minDistanceSegment = segment;
510 minDistance = minDistanceNew;
511 }
else if (minDistanceNew < minDistance) {
512 minDistance = minDistanceNew;
513 minDistanceSegment = segment;
522 strip = (minDistanceSegment - 1) * nStripsSegment;
524 [sector - 1][plane - 1][strip] * intersectionClhep;
525 y = intersectionLocal.y() / CLHEP::cm *
Unit::cm;
526 stripSegment = ceil((y - 0.5 * stripWidth) / stripWidth);
527 if (stripSegment <= 0)
529 else if (stripSegment > nStripsSegment)
530 stripSegment = nStripsSegment;
531 strip = stripSegment + (minDistanceSegment - 1) * nStripsSegment;
533 [sector - 1][plane - 1][strip - 1] * intersectionClhep;
534 x = intersectionLocal.x();
535 l =
m_GeoDat->getStripLength(strip);
540 if (fabs(x) > 0.5 * l)
543 section, layer, sector, plane, strip);
545 *strip1 = stripGlobal;
547 *strip2 = stripGlobal;