15 #include "GeoModelKernel/GeoDefinitions.h"
16 #include "GeoModelKernel/GeoTube.h"
20 #include "GeoModelKernel/throwExcept.h"
21 #include "GeoModelHelpers/TransformToStringConverter.h"
22 #include "GeoModelHelpers/GeoShapeUtils.h"
39 #define verbose_Bline false
71 std::vector<const Trk::Surface*> elementSurfaces;
75 if (
s) elementSurfaces.push_back(
s.get());
77 return elementSurfaces;
82 if (tubeLayer >
getNLayers() || tubeLayer < 1)
return false;
88 if (tubeLayer >
getNLayers() || tubeLayer < 1)
return false;
104 double nominalTubeLength = 0.;
109 if (istep < 0 || istep >=
m_nsteps) {
111 <<
" DEid = " <<
idHelperSvc()->toStringDetEl(
identify()) <<
" called with: tubeL, tube " << tubeLayer <<
" " <<
tube
118 double tlength = nominalTubeLength;
121 if (
manager()->MinimalGeoFlag() == 0) {
123 <<
" has cutouts, check for real tube length for tubeLayer, tube " << tubeLayer <<
" " <<
tube );
124 PVConstLink cv = getMaterialGeom();
125 int nGrandchildren = cv->getNChildVols();
126 if (nGrandchildren <= 0)
return tlength;
135 int packed_id =
tube + maxNTubesPerLayer * tubeLayer;
140 if (!
found &&
id == packed_id) {
147 ATH_MSG_DEBUG(
" MdtReadoutElement tube match found for BMG - input : tube(" <<
tube <<
"), layer("
148 << tubeLayer <<
") - output match : tube(" << ii % maxNTubesPerLayer <<
"), layer(" << ii / maxNTubesPerLayer
152 if (ii >= nGrandchildren) {
154 <<
" --- getTubeLength is looking for child with index ii=" << ii <<
" for tubeL and tube = " << tubeLayer <<
" "
161 <<
" --- getTubeLength is looking for child with index ii=" << ii <<
" for tubeL and tube = " << tubeLayer <<
" "
166 PVConstLink physChild = cv->getChildVol(ii);
167 const GeoShape* shape = physChild->getLogVol()->getShape();
168 if (shape ==
nullptr)
return tlength;
169 const GeoTube* theTube =
dynamic_cast<const GeoTube*
>(shape);
170 if (theTube !=
nullptr)
171 tlength = 2. * theTube->getZHalfLength();
174 <<
" out of (tubeLayer-1)*m_ntubesperlayer+tube with tl=" << tubeLayer <<
" tubes/lay=" <<
m_ntubesperlayer
177 if (std::abs(tlength - nominalTubeLength) > 0.1) {
179 <<
" is affected by a cutout: eff. length = " << tlength <<
" while nominal = " << nominalTubeLength
183 <<
" is NOT affected by the cutout: eff. length = " << tlength <<
" while nominal = " << nominalTubeLength);
196 double scalprod = c_ro.dot(x_ro);
228 int amdb_plus_minus1 = 1;
231 if (std::abs(
ms->xAmdbCRO()) > 10.) {
236 if (ROtubeFrame.z() < 0)
275 if (amdb_plus_minus1 == 0) {
276 ATH_MSG_WARNING(
"Unable to get the sign of RO side; signedRODistancefromTubeCenter returns 0" );
297 <<
"/" << tubeLayer <<
"/" <<
tube );
298 double xtube{0.}, ytube{0.}, ztube{0.};
310 PVConstLink cv = getMaterialGeom();
311 int nGrandchildren = cv->getNChildVols();
321 int packed_id =
tube + maxNTubesPerLayer * tubeLayer;
326 if (!
found &&
id == packed_id) {
334 ATH_MSG_DEBUG(
" MdtReadoutElement tube match found for BMG - input : tube(" <<
tube <<
"), layer("
335 << tubeLayer <<
") - output match : tube(" << ii % maxNTubesPerLayer <<
"), layer(" << ii / maxNTubesPerLayer
340 PVConstLink tv = cv->getChildVol(ii);
341 constexpr
double maxtol = 0.0000001;
343 if (std::abs(xtube - tubeTrans(0, 3)) > maxtol || std::abs(ztube - tubeTrans(2, 3)) > maxtol) {
344 THROW_EXCEPTION(__func__<<
"("<<tubeLayer<<
","<<
tube<<
")"<<
" - mismatch between local from tube-id/pitch/cutout position"<<
346 idHelperSvc()->toStringDetEl(
identify()) <<
"There are "<<nGrandchildren<<
" child volumes and "<<
350 if (tubeTrans(1, 3) > maxtol) {
353 <<
"/" << tubeLayer <<
"/" <<
tube);
356 THROW_EXCEPTION(__func__<<
"("<<tubeLayer<<
","<<
tube<<
")"<<
" - mismatch between local from tube-id/pitch/cutout position"<<
358 idHelperSvc()->toStringDetEl(
identify()) <<
"There are "<<nGrandchildren<<
" child volumes and "<<
412 if (istep < 0 || istep >=
m_nsteps) {
414 <<
" called with: tubeL, tube " << tubeLayer
422 double xtube{0.}, ztube{0.};
436 ATH_MSG_WARNING(__func__<<
"() :"<<__LINE__<<
" called with tubeLayer or tube out of range in chamber "
438 <<
" max " <<
m_ntubesperlayer <<
" will compute deformation for first tube in this chamber" );
446 ptr.set(std::make_unique<Amg::Transform3D>(std::move(trans)));
464 if (!
ms->hasBLines() && !
ms->hasMdtAsBuiltParams()) {
465 return Amg::Transform3D::Identity();
472 double height =
barrel() ?
ms->ZsizeMdtStation() :
ms->RsizeMdtStation();
473 double thickness =
barrel() ?
ms->RsizeMdtStation() :
ms->ZsizeMdtStation();
476 <<
", layer: "<<tubeLayer<<
", tube: "<<
tube<<
", fixedPoint: "<<
Amg::toString(fixedPoint)<<
" / "
477 <<
Amg::toString(
ms->getGeoTransform()->getDefTransform() *
ms->getNativeToAmdbLRS().inverse()* fixedPoint)
478 <<
", height: "<<height<<
", thickness: "<<thickness
484 if (std::abs(height - heightML) > 10.) {
486 <<
"Different ML and MDTStation length in the precision coord.direction --- MultiLayerHeight, MDTstationHeigh "
487 << heightML <<
" " << height <<
" MultiLayerThickness, MDTstationThickness " << thicknessML <<
" " << thickness);
503 Amg::Vector3D pt_end1 = toAMDB * pt_center + halftubelen * Amg::Vector3D::UnitX();
504 Amg::Vector3D pt_end2 = toAMDB * pt_center - halftubelen * Amg::Vector3D::UnitX();
510 if (
ms->hasMdtAsBuiltParams()) {
515 if (
ms->hasBLines()) {
519 pt_end1_new =
posOnDefChamWire(pt_end1_new, moduleWidthS, moduleWidthL, height, thickness, fixedPoint);
522 pt_end2_new =
posOnDefChamWire(pt_end2_new, moduleWidthS, moduleWidthL, height, thickness, fixedPoint);
526 pt_end1 = fromAMDB * pt_end1;
527 pt_end2 = fromAMDB * pt_end2;
528 pt_end1_new = fromAMDB * pt_end1_new;
529 pt_end2_new = fromAMDB * pt_end2_new;
536 const Amg::Vector3D pt_center_new = 0.5 * (pt_end1_new + pt_end2_new);
541 const Amg::Vector3D rotation_vector = old_direction.cross(new_direction);
545 if (rotation_vector.mag() > 10. * std::numeric_limits<double>::epsilon()) {
546 const Amg::AngleAxis3D wire_rotation(std::asin(rotation_vector.mag()), rotation_vector.unit());
579 const double moduleWidthS,
580 const double moduleWidthL,
582 const double thickness,
622 ATH_MSG_VERBOSE(
"** In "<<__func__<<
" - moduleWidthS " << moduleWidthS<<
", moduleWidthL: "
623 <<moduleWidthL <<
", height: "<<height<<
", thickness: " <<thickness <<
"." );
625 <<
" in the amdb-szt frame" );
627 double s0mdt = locAMDBPos.x();
628 if (std::abs(fixedPoint.x()) > 0.01) {
629 s0mdt = locAMDBPos.x() - fixedPoint.x();
631 double z0mdt = locAMDBPos.y();
633 if (std::abs(fixedPoint.y()) > 0.01) {
634 z0mdt = locAMDBPos.y() - fixedPoint.y();
636 double t0mdt = locAMDBPos.z();
638 if (std::abs(fixedPoint.z()) > 0.01) t0mdt = locAMDBPos.z() - fixedPoint.z();
639 if (z0mdt < 0 || t0mdt < 0) {
640 ATH_MSG_WARNING(
""<<__func__<<
": correcting the local position of a point outside the mdt station (2 multilayers) volume -- RE "
644 ATH_MSG_VERBOSE(
"** In "<<__func__<<
" - correct for offset of B-line fixed point " << s0mdt <<
" " << z0mdt <<
" " << t0mdt);
646 double ds{0.},dz{0.},
dt{0.};
647 double width_actual = moduleWidthS + (moduleWidthL - moduleWidthS) * (z0mdt / height);
648 double s_rel = s0mdt / (width_actual / 2.);
649 double z_rel = (z0mdt - height / 2.) / (height / 2.);
650 double t_rel = (t0mdt - thickness / 2.) / (thickness / 2.);
652 ATH_MSG_VERBOSE(
"** In "<<__func__<<
" - width_actual: "<<width_actual<<
", s_rel: "<<s_rel<<
", z_rel: "
653 <<z_rel<<
", t_rel: "<<t_rel);
656 if ((sp != 0) || (sn != 0)) {
657 double ztmp = z_rel * z_rel - 1;
658 dt += 0.5 * (sp + sn) * ztmp + 0.5 * (sp - sn) * ztmp * s_rel;
663 dt -= tw * s_rel * z_rel;
664 dz += tw * s_rel * t_rel * thickness / height;
665 ATH_MSG_VERBOSE(
"** In "<<__func__<<
": tw=" << tw <<
" dt, dz " <<
dt <<
" " << dz );
670 double egppm =
eg / 1000.;
682 if ((ep != 0) || (
en != 0)) {
685 double phi = 0.5 * (ep +
en) * s_rel * s_rel + 0.5 * (ep -
en) * s_rel;
686 double localDt = phi * (t0mdt - thickness / 2.);
687 double localDz = phi * (z0mdt - height / 2.);
693 deformPos[0] = locAMDBPos[0] +
ds;
694 deformPos[1] = locAMDBPos[1] + dz;
695 deformPos[2] = locAMDBPos[2] +
dt;
713 const int tubeLayer,
const int tube)
const {
721 <<
" tubeLayer=" << tubeLayer <<
" tube=" <<
tube);
723 constexpr
int nsid = 2;
726 std::array<MdtAsBuiltPar::tubeSide_t, nsid> tubeSide{tubeSide_t::POS, tubeSide_t::NEG};
727 std::array<Amg::Vector3D, nsid> wireEnd{locAMDBWireEndP, locAMDBWireEndN};
728 multilayer_t ml = (
getMultilayer() == 1) ? multilayer_t::ML1 : multilayer_t::ML2;
731 const int ref_layer = (ml == multilayer_t::ML1) ?
m_nlayers : 1;
734 for (
int isid = 0; isid < nsid; ++isid) {
736 double xref{0.}, yref{0.}, zref{0.};
747 reference_point = toAMDB * reference_point;
757 int layer_delta = tubeLayer;
758 if (ml == multilayer_t::ML1) layer_delta =
m_nlayers + 1 - tubeLayer;
761 double zpitch =
params->zpitch(ml, tubeSide[isid]);
762 double ypitch =
params->ypitch(ml, tubeSide[isid]);
763 double stagg = (
double)
params->stagg(ml, tubeSide[isid]);
765 double y0 =
params->y0(ml, tubeSide[isid]);
766 double z0 =
params->z0(ml, tubeSide[isid]);
769 int do_stagg = (layer_delta - 1) % 2;
770 double offset_stagg = ((
double)do_stagg) * 0.5 * zpitch * stagg;
771 Amg::Vector3D end_plug(0., (
tube - 1.0) * zpitch + offset_stagg, (layer_delta - 1) * ypitch);
774 end_plug = amgTransf * end_plug;
783 Amg::Vector3D ret(reference_point.x() + xshift, reference_point.y() +
z0 + end_plug.y(),
784 reference_point.z() + y0 + end_plug.z());
790 << isid <<
", Delta " <<
Amg::toString(ret - wireEnd[isid]) );
792 ATH_MSG_VERBOSE(( tubeSide[isid] == tubeSide_t::POS ?
"positive" :
"negative")<<
" wire end has moved from "
796 if (tubeSide[isid] == tubeSide_t::POS)
797 locAMDBWireEndP = ret;
799 locAMDBWireEndN = ret;
804 return std::make_unique<GeoInfo>(std::move(
transform));
810 ATH_MSG_WARNING(__func__<<
"() :"<<__LINE__<<
" called with tubeLayer or tube out of range in chamber "
812 <<
" max " <<
m_ntubesperlayer <<
" will compute transform for first tube in this chamber" );
834 if (itube >= ntot_tubes) {
835 ATH_MSG_WARNING(
"surface called with tubeLayer or tube out of range in chamber "
837 <<
" max " <<
m_ntubesperlayer <<
" will compute surface for first tube in this chamber" );
845 ptr.set(std::make_unique<Trk::StraightLineSurface>(*
this,
id));
860 if (istep < 0 || istep >= ntot_steps) {
862 <<
" DEid = " <<
idHelperSvc()->toStringDetEl(
identify()) <<
" called with: tubeL, tube " << tubeLayer <<
" " <<
tube
864 ATH_MSG_WARNING(
"Please run in DEBUG mode to get extra diagnostic; setting istep = 0" );
889 surfaceTRotation.col(0) = muonTRotation.col(1);
890 surfaceTRotation.col(1) = muonTRotation.col(2);
891 surfaceTRotation.col(2) = muonTRotation.col(0);
894 trans3D.pretranslate(
transform().translation());
970 ATH_MSG_VERBOSE(
"global Surface / Bounds pointers " << tmpSurface <<
" " << tmpBounds );
983 PVConstLink cv = getMaterialGeom();
987 int packed_id =
tube + maxNTubesPerLayer *
tl;
991 if (!
found &&
id == packed_id) {
1026 if (layer < 1 || layer >
getNLayers())
return false;