ATLAS Offline Software
Loading...
Searching...
No Matches
MuonDetDescr/MuonReadoutGeometry/src/sTgcReadoutElement.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5/***************************************************************************
6 The sTgc detector = an assembly module = STGC in amdb
7 ----------------------------------------------------
8***************************************************************************/
9
11
14
15#include <GeoModelKernel/GeoLogVol.h>
16#include <GeoModelKernel/GeoDefinitions.h>
17#include <GeoModelHelpers/StringUtils.h>
18#include <GeoModelHelpers/TransformToStringConverter.h>
19#include <GeoModelHelpers/getChildNodesWithTrf.h>
20
21#include <cmath>
22#include <ext/alloc_traits.h>
23#include <map>
24#include <memory>
25#include <stdexcept>
26#include <utility>
27
28#include "GeoModelKernel/GeoFullPhysVol.h"
29#include "GaudiKernel/SystemOfUnits.h"
41
42#include "GaudiKernel/ISvcLocator.h"
48
49#define THROW_EXCEPTION_RE(MSG) \
50 { \
51 std::stringstream sstr{}; \
52 sstr<<"sTgcReadoutElement - "<<idHelperSvc()->toStringDetEl(identify())<<" "<<__LINE__<<": "; \
53 sstr<<MSG; \
54 throw std::runtime_error(sstr.str()); \
55 } \
56
57using namespace GeoStrUtils;
58namespace MuonGM {
59
60 //============================================================================
61 sTgcReadoutElement::sTgcReadoutElement(GeoVFullPhysVol* pv, const std::string& stName, int zi, int fi, int mL, MuonDetectorManager* mgr)
62 : MuonClusterReadoutElement(pv, mgr, Trk::DetectorElemType::sTgc)
63 , m_ml(mL) {
64
65 std::string fixName = (stName[1] == 'L') ? "STL" : "STS";
66 Identifier id = mgr->stgcIdHelper()->channelID(fixName, zi, fi, mL, 1,
68
69 setStationName(fixName);
71 setIdentifier(id); // representative identifier, with stName, stEta, stPhi, mL
72 }
73
74
75 //============================================================================
77
78 //============================================================================
80
81 PVConstLink pvc {getMaterialGeom()};
82 auto sensitiveVol = getAllSubVolumes(pvc,[](const GeoChildNodeWithTrf& node){
83 return node.nodeName.find("Gas") != std::string::npos;
84 });
85 assert(sensitiveVol.size() == m_nlayers);
86 for (unsigned int llay = 0; llay< sensitiveVol.size(); ++llay) {
87 m_Xlg[llay] = sensitiveVol[llay].transform;
88 }
89
90 SmartIF<IGeoDbTagSvc> geoDbTag{Gaudi::svcLocator()->service("GeoDbTagSvc")};
91 SmartIF<IRDBAccessSvc> accessSvc{Gaudi::svcLocator()->service(geoDbTag->getParamSvcName())};
92
93 IRDBRecordset_ptr nswdimRec = accessSvc->getRecordsetPtr("NSWDIM","","");
94 IRDBRecordset_ptr wstgcRec = accessSvc->getRecordsetPtr("WSTGC","","");
95 IRDBRecordset_ptr nswPars = accessSvc->getRecordsetPtr("NSWPARS","","");
96
97 PVConstLink parent = getMaterialGeom()->getParent();
98 unsigned int index=parent->indexOf(getMaterialGeom()).value();
99 std::string pVName=parent->getNameOfChildVol(index);
100 float yCutoutCathode(0);
101 if (nswPars->size()==0) {
102 THROW_EXCEPTION_RE("Error, cannot access NSWPARS record!");
103 } else {
104 yCutoutCathode=(*nswPars)[0]->getFloat("NSW_sTGC_yCutoutCathode");
105 }
106
107 for (size_t w=0;w<nswdimRec->size();w++) {
108 const IRDBRecord *nswdim = (*nswdimRec)[w];
109 const std::string type = nswdim->getString("NSW_TYPE").substr(5,4);
110 std::string logVolSubName=getMaterialGeom()->getLogVol()->getName().substr(7,4);
111 if (type==logVolSubName) {
112 setSsize(nswdim->getDouble("BASE_WIDTH")); // bottom base length (full chamber)
113 setLongSsize(nswdim->getDouble("TOP_WIDTH")); // top base length (full chamber)
114 setRsize(nswdim->getDouble("LENGTH")); // height of the trapezoid (full chamber)
115 break;
116 }
117 }
118
119 for (unsigned int ind = 0; ind < wstgcRec->size(); ind++) {
120 std::string WSTGC_TYPE = (*wstgcRec)[ind]->getString("WSTGC_TYPE");
121 if (WSTGC_TYPE.size()<=8)
122 THROW_EXCEPTION_RE("Malformed WSTGC_TYPE = " << WSTGC_TYPE);
123 if (getStationName()[2] != WSTGC_TYPE[6]) continue;
124 if (std::abs(getStationEta()) != static_cast<int>(WSTGC_TYPE[7]-'0')) continue;
125 if (getStationName()[2] == 'S' && WSTGC_TYPE[8] != (m_ml ==2 ?'P' : 'C')) continue;
126 if (getStationName()[2] == 'L' && WSTGC_TYPE[8] != (m_ml ==2 ?'C' : 'P')) continue;
127
128 const double gasTck = (*wstgcRec)[ind]->getDouble("gasTck");
129 const double Tck = (*wstgcRec)[ind]->getDouble("Tck");
130 const double xFrame = (*wstgcRec)[ind]->getDouble("xFrame");
131 const double ylFrame = (*wstgcRec)[ind]->getDouble("ylFrame");
132 const double ysFrame = (*wstgcRec)[ind]->getDouble("ysFrame");
133 const double wirePitch = (*wstgcRec)[ind]->getDouble("wirePitch");
134 const double stripPitch = (*wstgcRec)[ind]->getDouble("stripPitch");
135 const double stripWidth = (*wstgcRec)[ind]->getDouble("stripWidth");
136 const double sPadWidth = (*wstgcRec)[ind]->getDouble("sPadWidth");
137 const double lPadWidth = (*wstgcRec)[ind]->getDouble("lPadWidth");
138 const double anglePadPhi = (*wstgcRec)[ind]->getDouble("anglePadPhi");
139 const double sStripWidth = (*wstgcRec)[ind]->getDouble("sStripWidth");
140 const double lStripWidth = (*wstgcRec)[ind]->getDouble("lStripWidth");
141 const int wireGroupWidth = (*wstgcRec)[ind]->getInt("wireGroupWidth");
142 const int nStrips = (*wstgcRec)[ind]->getInt("nStrips");
143 const std::vector<double> padH = tokenizeDouble((*wstgcRec)[ind]->getString("padH"),";");
144 const std::vector<int> nPadPhi = tokenizeInt((*wstgcRec)[ind]->getString("nPadPhi"),";");
145 const std::vector<double> firstPadPhiDivision_A = tokenizeDouble((*wstgcRec)[ind]->getString("firstPadPhiDivision_A"),";");
146 const std::vector<double> PadPhiShift_A = tokenizeDouble((*wstgcRec)[ind]->getString("PadPhiShift_A"),";");
147 const std::vector<int> nPadH = tokenizeInt((*wstgcRec)[ind]->getString("nPadH"),";");
148 const std::vector<double> firstPadH = tokenizeDouble((*wstgcRec)[ind]->getString("firstPadH"),";");
149 const std::vector<double> firstPadRow = tokenizeDouble((*wstgcRec)[ind]->getString("firstPadRow"),";");
150 const std::vector<double> wireCutout = tokenizeDouble((*wstgcRec)[ind]->getString("wireCutout"),";");
151 const std::vector<int> nWires = tokenizeInt((*wstgcRec)[ind]->getString("nWires"),";");
152 const std::vector<int> firstWire = tokenizeInt((*wstgcRec)[ind]->getString("firstWire"),";");
153 const std::vector<double> firstStripWidth = tokenizeDouble((*wstgcRec)[ind]->getString("firstStripWidth"),";");
154 const std::vector<int> nWireGroups = tokenizeInt((*wstgcRec)[ind]->getString("nWireGroups"),";");
155 const std::vector<double> firstWireGroup = tokenizeDouble((*wstgcRec)[ind]->getString("firstWireGroup"),";");
156
157 char sector_l = getStationName()[2];
158 int stEta = std::abs(getStationEta());
159 int Etasign = getStationEta() / stEta;
160 std::string side = (Etasign > 0) ? "A" : "C";
161 m_diamondShape = sector_l == 'L' && stEta == 3;
162
163
164 // Get frame widths
165 setZsize(Tck); // thickness (full chamber)
166
167 double yCutout = m_diamondShape ? yCutoutCathode: 0.0; // y of cutout of trapezoid (only in outermost detectors)
168
169
170 // Radial shift of the local frame origin w.r.t. the center of the quadruplet.
171 // For diamond shape (QL3) the origin is on the cutout base. For the rest, the it is at the center
172 // of the active area, therefore the shift is half the difference of the top and bottom frame widths.
173 m_offset = (m_diamondShape) ? 0.5*getRsize() - (yCutout + ylFrame) : -0.5*(ylFrame - ysFrame);
174
175 //-------------------
176 // Strips
177 //-------------------
178 for (int il = 0; il < m_nlayers; il++) {
181 if (yCutout == 0.) {
182 m_etaDesign[il].defineTrapezoid(0.5 * sStripWidth,
183 0.5 * lStripWidth,
184 0.5 * (getRsize() - ysFrame - ylFrame));
185 } else {
186 m_etaDesign[il].defineDiamond(0.5 * sStripWidth,
187 0.5 * lStripWidth,
188 0.5 * (getRsize() - ysFrame - ylFrame), yCutout);
189 }
190 m_etaDesign[il].inputPitch = stripPitch;
191 m_etaDesign[il].inputWidth = stripWidth;
192 m_etaDesign[il].thickness = gasTck;
193 m_etaDesign[il].firstPitch = firstStripWidth[il];
194 m_etaDesign[il].setFirstPos(m_diamondShape ? -(m_etaDesign[il].xSize()- yCutout) + m_etaDesign[il].firstPitch
195 : -0.5 * m_etaDesign[il].xSize()+ m_etaDesign[il].firstPitch);
196 m_etaDesign[il].nch = nStrips;
197 }
198
199 //-------------------
200 // Wires
201 //-------------------
202 for (int il = 0; il < m_nlayers; il++) {
205 if (yCutout == 0.) {
206 m_phiDesign[il].defineTrapezoid(0.5 * sPadWidth,
207 0.5 * lPadWidth,
208 0.5 * (getRsize() - ysFrame - ylFrame) );
209 } else {
210 m_phiDesign[il].defineDiamond(0.5 * sPadWidth,
211 0.5 * lPadWidth,
212 0.5 * (getRsize() - ysFrame - ylFrame), yCutout);
213 }
214 m_phiDesign[il].inputPitch = wirePitch;
215 m_phiDesign[il].inputWidth = 0.015;
216 m_phiDesign[il].thickness = getZsize();
217 m_phiDesign[il].setFirstPos(firstWire[il]); // Position of 1st wire, accounts for staggering
218 m_phiDesign[il].firstPitch = firstWireGroup[il]; // Number of Wires in 1st group, group staggering
219 m_phiDesign[il].groupWidth = wireGroupWidth; // Number of Wires normal group
220 m_phiDesign[il].nGroups = nWireGroups[il]; // Number of Wire Groups
221 m_phiDesign[il].wireCutout = wireCutout[il]; // Size of "active" wire region for digits
222 m_phiDesign[il].nch = nWires[il];
223
224 }
225
226 //-------------------
227 // Pads
228 //-------------------
229 double radius = absTransform().translation().perp() + m_offset;
230 for (int il = 0; il < m_nlayers; il++) {
231 m_padDesign[il].Length = getRsize();
232 m_padDesign[il].sWidth = getSsize();
233 m_padDesign[il].lWidth = getLongSsize();
234 m_padDesign[il].Size = getRsize() - ylFrame - ysFrame;
235 m_padDesign[il].xFrame = xFrame;
236 m_padDesign[il].ysFrame = ysFrame;
237 m_padDesign[il].ylFrame = ylFrame;
238 m_padDesign[il].yCutout = yCutout;
239 m_padDesign[il].etasign = Etasign;
240 m_padDesign[il].setR(radius);
241 m_padDesign[il].sPadWidth = sPadWidth;
242 m_padDesign[il].lPadWidth = lPadWidth;
243
244 m_padDesign[il].nPadColumns = nPadPhi[il];
245
246 // The C side of the NSW is mirrored instead of rotated
247 // We should be using the same values for the pads for both A and C
248 // It is easier for us to simply read the same correct value once
249 // whereas changing the XML and the reading functions will make this incompatible with past versions
250 // Alexandre Laurier 12 Sept 2018
251 m_padDesign[il].firstPhiPos= firstPadPhiDivision_A[il];
252 m_padDesign[il].inputPhiPitch = anglePadPhi; // stEta<2 ? PAD_PHI_DIVISION/PAD_PHI_SUBDIVISION : PAD_PHI_DIVISION ;
253 m_padDesign[il].PadPhiShift= PadPhiShift_A[il];
254 m_padDesign[il].padEtaMin = firstPadRow[il]; // FIRST_PAD_ROW_DIVISION[2*sector+(m_ml-1)][stEta-1][il];
255 m_padDesign[il].nPadH = nPadH[il];
256 m_padDesign[il].padEtaMax = m_padDesign[il].padEtaMin + m_padDesign[il].nPadH; // PAD_ROWS[2*sector+(m_ml-1)][stEta-1][il];
257 m_padDesign[il].firstRowPos = firstPadH[il]; // H_PAD_ROW_0[2*sector+(m_ml-1)][il];
258 m_padDesign[il].inputRowPitch = padH[il]; // PAD_HEIGHT[2*sector+(m_ml-1)][il];
259
260 if (sector_l == 'L') {
261 m_padDesign[il].isLargeSector = 1;
262 m_padDesign[il].sectorOpeningAngle = m_padDesign[il].largeSectorOpeningAngle;
263 } else {
264 m_padDesign[il].isLargeSector = 0;
265 m_padDesign[il].sectorOpeningAngle = m_padDesign[il].smallSectorOpeningAngle;
266 }
267 m_padDesign[il].thickness = thickness;
268 ATH_MSG_DEBUG("initDesign: " << idHelperSvc()->toStringDetEl(identify()) << " layer " << il << ", pad phi angular width "
269 << m_padDesign[il].inputPhiPitch << ", eta pad size " << m_padDesign[il].inputRowPitch
270 << " Length: " << m_padDesign[il].Length << " sWidth: " << m_padDesign[il].sWidth
271 << " lWidth: " << m_padDesign[il].lWidth << " firstPhiPos:" << m_padDesign[il].firstPhiPos
272 << " padEtaMin:" << m_padDesign[il].padEtaMin << " padEtaMax:" << m_padDesign[il].padEtaMax
273 << " firstRowPos:" << m_padDesign[il].firstRowPos << " inputRowPitch:" << m_padDesign[il].inputRowPitch
274 << " thickness:" << m_padDesign[il].thickness << " sPadWidth: " << m_padDesign[il].sPadWidth
275 << " lPadWidth: " << m_padDesign[il].lPadWidth << " xFrame: " << m_padDesign[il].xFrame
276 << " ysFrame: " << m_padDesign[il].ysFrame << " ylFrame: " << m_padDesign[il].ylFrame
277 << " yCutout: " << m_padDesign[il].yCutout );
278 }
279 }
280 }
282
283 if (manager()->MinimalGeoFlag() == 0) {
284 PVConstLink pvc {getMaterialGeom()};
285 unsigned int nchildvol = pvc->getNChildVols();
286 int llay = 0;
287 std::string::size_type npos;
288 for (unsigned ich = 0; ich < nchildvol; ++ich) {
289 PVConstLink pc = pvc->getChildVol(ich);
290 std::string childname = (pc->getLogVol())->getName();
291
292 ATH_MSG_DEBUG("Volume Type: " << pc->getLogVol()->getShape()->type());
293 if ((npos = childname.find("Sensitive")) == std::string::npos) {
294 continue;
295 }
296 ++llay;
297 if (llay > 4) {
298 ATH_MSG_DEBUG("number of sTGC layers > 4: increase transform array size");
299 continue;
300 }
301 m_Xlg[llay - 1] = pvc->getXToChildVol(ich);
302 }
303 assert(m_nlayers == llay);
304 }
305 char sector_l = getStationName().substr(2, 1) == "L" ? 'L' : 'S';
306 int stEta = std::abs(getStationEta());
307 int Etasign = getStationEta() / stEta;
308 std::string side = (Etasign > 0) ? "A" : "C";
309 m_diamondShape = sector_l == 'L' && stEta == 3;
310
311 ATH_MSG_DEBUG("station name" << getStationName());
312
313 sTGCDetectorHelper aHelper;
314 sTGCDetectorDescription* stgc = aHelper.Get_sTGCDetector(sector_l, stEta, getStationPhi(), m_ml, side.back());
315
316 ATH_MSG_DEBUG( "Found sTGC Detector " << stgc->GetName() );
317
319 if (!tech) THROW_EXCEPTION_RE(" Failed To get Technology for stgc element:"<< stgc->GetName());
320
321
322 // Get Chamber length, width and frame widths
323 setSsize(stgc->sWidth()); // bottom base length (full chamber)
324 setLongSsize(stgc->lWidth()); // top base length (full chamber)
325 setRsize(stgc->Length()); // height of the trapezoid (full chamber)
326 setZsize(stgc->Tck()); // thickness (full chamber)
327 double ysFrame = stgc->ysFrame(); // Frame thickness on short parallel edge
328 double ylFrame = stgc->ylFrame(); // Frame thickness on long parallel edge
329 double xFrame = stgc->xFrame(); // Frame thickness of non parallel edges
330 double yCutout = stgc->yCutoutCathode(); // y of cutout of trapezoid (only in outermost detectors)
332
333
334 // Radial shift of the local frame origin w.r.t. the center of the quadruplet.
335 // For diamond shape (QL3) the origin is on the cutout base. For the rest, the it is at the center
336 // of the active area, therefore the shift is half the difference of the top and bottom frame widths.
337 m_offset = (m_diamondShape) ? 0.5*getRsize() - (yCutout + ylFrame) : -0.5*(ylFrame - ysFrame);
338
339 //-------------------
340 // Strips
341 //-------------------
342 for (int il = 0; il < m_nlayers; il++) {
343 // identifier of the first channel - strip plane - to retrieve max number of strips
346 if (yCutout == 0.) {
347 m_etaDesign[il].defineTrapezoid(0.5 * roParam.sStripWidth,
348 0.5 * roParam.lStripWidth,
349 0.5 * (getRsize() - ysFrame - ylFrame));
350 } else {
351 m_etaDesign[il].defineDiamond(0.5 * roParam.sStripWidth,
352 0.5 * roParam.lStripWidth,
353 0.5 * (getRsize() - ysFrame - ylFrame), yCutout);
354 }
355 m_etaDesign[il].inputPitch = stgc->stripPitch();
356 m_etaDesign[il].inputWidth = stgc->stripWidth();
357 m_etaDesign[il].thickness = tech->gasThickness;
358 m_etaDesign[il].firstPitch = roParam.firstStripWidth[il];
359 m_etaDesign[il].setFirstPos((m_diamondShape) ? -(m_etaDesign[il].xSize()- yCutout) + m_etaDesign[il].firstPitch
360 : -0.5 * m_etaDesign[il].xSize()+ m_etaDesign[il].firstPitch);
361 m_etaDesign[il].nch = roParam.nStrips;
362
363 ATH_MSG_DEBUG("initDesign:" << getStationName() << " layer " << il
364 << ", strip pitch " << m_etaDesign[il].inputPitch
365 << ", nstrips " << m_etaDesign[il].nch
366 << ", firstPos: " << m_etaDesign[il].firstPos() );
367 }
368
369 //-------------------
370 // Wires
371 //-------------------
372 for (int il = 0; il < m_nlayers; il++) {
375 if (yCutout == 0.) {
376 m_phiDesign[il].defineTrapezoid(0.5 * roParam.sPadWidth,
377 0.5 * roParam.lPadWidth,
378 0.5 * (getRsize() - ysFrame - ylFrame) );
379 } else {
380 m_phiDesign[il].defineDiamond(0.5 * roParam.sPadWidth,
381 0.5 * roParam.lPadWidth,
382 0.5 * (getRsize() - ysFrame - ylFrame), yCutout);
383 }
384 m_phiDesign[il].inputPitch = stgc->wirePitch();
385 m_phiDesign[il].inputWidth = 0.015;
386 m_phiDesign[il].thickness = getZsize();
387 m_phiDesign[il].setFirstPos(roParam.firstWire[il]); // Position of 1st wire, accounts for staggering
388 m_phiDesign[il].firstPitch = roParam.firstWireGroup[il]; // Number of Wires in 1st group, group staggering
389 m_phiDesign[il].groupWidth = roParam.wireGroupWidth; // Number of Wires normal group
390 m_phiDesign[il].nGroups = roParam.nWireGroups[il]; // Number of Wire Groups
391 m_phiDesign[il].wireCutout = roParam.wireCutout[il]; // Size of "active" wire region for digits
392 m_phiDesign[il].nch = roParam.nWires[il];
393
394 ATH_MSG_DEBUG( "initDesign:" << getStationName() << " layer " << il << ", wireGang pitch "
395 << m_phiDesign[il].inputPitch << ", nWireGangs " << m_phiDesign[il].nch );
396 }
397
398 //-------------------
399 // Pads
400 //-------------------
401 double radius = absTransform().translation().perp() + m_offset;
402 for (int il = 0; il < m_nlayers; il++) {
403 m_padDesign[il].Length = getRsize();
404 m_padDesign[il].sWidth = getSsize();
405 m_padDesign[il].lWidth = getLongSsize();
406 m_padDesign[il].Size = getRsize() - ylFrame - ysFrame;
407 m_padDesign[il].xFrame = xFrame;
408 m_padDesign[il].ysFrame = ysFrame;
409 m_padDesign[il].ylFrame = ylFrame;
410 m_padDesign[il].yCutout = yCutout;
411 m_padDesign[il].etasign = Etasign;
412 m_padDesign[il].setR(radius);
413 m_padDesign[il].sPadWidth = roParam.sPadWidth;
414 m_padDesign[il].lPadWidth = roParam.lPadWidth;
415 m_padDesign[il].nPadColumns = roParam.nPadPhi[il];
416
417 // The C side of the NSW is mirrored instead of rotated
418 // We should be using the same values for the pads for both A and C
419 // It is easier for us to simply read the same correct value once
420 // whereas changing the XML and the reading functions will make this incompatible with past versions
421 // Alexandre Laurier 12 Sept 2018
422 m_padDesign[il].firstPhiPos = roParam.firstPadPhiDivision_A[il];
423 m_padDesign[il].inputPhiPitch = roParam.anglePadPhi; // stEta<2 ? PAD_PHI_DIVISION/PAD_PHI_SUBDIVISION : PAD_PHI_DIVISION ;
424 m_padDesign[il].PadPhiShift = roParam.PadPhiShift_A[il];
425 m_padDesign[il].padEtaMin = roParam.firstPadRow[il]; // FIRST_PAD_ROW_DIVISION[2*sector+(m_ml-1)][stEta-1][il];
426 m_padDesign[il].nPadH = roParam.nPadH[il];
427 m_padDesign[il].padEtaMax = m_padDesign[il].padEtaMin + roParam.nPadH[il]; // PAD_ROWS[2*sector+(m_ml-1)][stEta-1][il];
428 m_padDesign[il].firstRowPos = roParam.firstPadH[il]; // H_PAD_ROW_0[2*sector+(m_ml-1)][il];
429 m_padDesign[il].inputRowPitch = roParam.padH[il]; // PAD_HEIGHT[2*sector+(m_ml-1)][il];
430
431 if (sector_l == 'L') {
432 m_padDesign[il].isLargeSector = 1;
433 m_padDesign[il].sectorOpeningAngle = m_padDesign[il].largeSectorOpeningAngle;
434 } else {
435 m_padDesign[il].isLargeSector = 0;
436 m_padDesign[il].sectorOpeningAngle = m_padDesign[il].smallSectorOpeningAngle;
437 }
438
439 m_padDesign[il].thickness = thickness;
440
441 ATH_MSG_DEBUG( "initDesign: " << idHelperSvc()->toStringDetEl(identify())
442 << " layer " << il<< ", pad phi angular width "
443 << m_padDesign[il].inputPhiPitch << ", eta pad size " << m_padDesign[il].inputRowPitch
444 << " Length: " << m_padDesign[il].Length << " sWidth: " << m_padDesign[il].sWidth
445 << " lWidth: " << m_padDesign[il].lWidth << " firstPhiPos:" << m_padDesign[il].firstPhiPos
446 << " padEtaMin:" << m_padDesign[il].padEtaMin << " padEtaMax:" << m_padDesign[il].padEtaMax
447 << " firstRowPos:" << m_padDesign[il].firstRowPos << " inputRowPitch:" << m_padDesign[il].inputRowPitch
448 << " thickness:" << m_padDesign[il].thickness << " sPadWidth: " << m_padDesign[il].sPadWidth
449 << " lPadWidth: " << m_padDesign[il].lPadWidth << " xFrame: " << m_padDesign[il].xFrame
450 << " ysFrame: " << m_padDesign[il].ysFrame << " ylFrame: " << m_padDesign[il].ylFrame
451 << " yCutout: " << m_padDesign[il].yCutout );
452 }
453 }
454
455 void sTgcReadoutElement::initDesign(double thickness) {
456
457
458
459 SmartIF<IGeoDbTagSvc> geoDbTag{Gaudi::svcLocator()->service("GeoDbTagSvc")};
460 if (!geoDbTag) THROW_EXCEPTION_RE( "Could not locate GeoDbTagSvc" );
461 if (geoDbTag->getSqliteReader()) initDesignFromSQLite(thickness);
462 else initDesignFromAGDD(thickness);
463
464 }
465
466 //============================================================================
468
469 if (m_surfaceData) {
470 ATH_MSG_WARNING("calling fillCache on an already filled cache" );
471 return;
472 }
473
474 m_surfaceData = std::make_unique<SurfaceData>();
475
476 for (int layer{0}; layer < m_nlayers; ++layer) {
477
478 // Define the geometry for the strips, pads and wires of this readout element.
479 // For QL3 (cutoff trapezoid), diamondBounds are used, while trapezoid bounds are used for the rest.
480 // The assigned coordinate along the layer normal is at the center of the gas gap;
481 // wires are considered at x=0, while strips and pads are shifted by +10/-10 microns.
482
483 //-------------------
484 // Layer boundaries
485 //-------------------
486
487 if (m_diamondShape) {
488 m_surfaceData->m_surfBounds.push_back(std::make_unique<Trk::RotatedDiamondBounds>(m_etaDesign[layer].minYSize() / 2.,
489 m_etaDesign[layer].maxYSize() / 2.,
490 m_etaDesign[layer].maxYSize() / 2.,
491 m_etaDesign[layer].xSize() / 2. - m_etaDesign[layer].yCutout() / 2,
492 m_etaDesign[layer].yCutout() / 2)); // strips
493
494 m_surfaceData->m_surfBounds.push_back(std::make_unique<Trk::DiamondBounds>(m_padDesign[layer].sPadWidth / 2.,
495 m_padDesign[layer].lPadWidth / 2.,
496 m_padDesign[layer].lPadWidth / 2.,
497 m_padDesign[layer].Size / 2. - m_padDesign[layer].yCutout / 2, m_padDesign[layer].yCutout / 2)); // pad and wires
498
499 } else {
500 m_surfaceData->m_surfBounds.push_back(std::make_unique<Trk::RotatedTrapezoidBounds>(m_etaDesign[layer].xSize() / 2.,
501 m_etaDesign[layer].minYSize() / 2.,
502 m_etaDesign[layer].maxYSize() / 2.)); // strips
503
504 m_surfaceData->m_surfBounds.push_back(std::make_unique<Trk::TrapezoidBounds>(m_padDesign[layer].sPadWidth /2.,
505 m_padDesign[layer].lPadWidth / 2.,
506 m_padDesign[layer].Size / 2.));
507
508
509 }
510
511 //-------------------
512 // Wires
513 //-------------------
514
515 // identifier of the first channel - wire plane - locX along phi, locY max->min R
516 Identifier id = m_idHelper.channelID(getStationName(), getStationEta(), getStationPhi(), m_ml, layer + 1, 2, 1);
517
518 m_surfaceData->m_layerSurfaces.push_back(std::make_unique<Trk::PlaneSurface>(*this, id));
519 m_surfaceData->m_layerTransforms.push_back(
520 absTransform() // transformation from chamber to ATLAS frame
521 * m_delta // transformations from the alignment group
522 * m_Xlg[layer] // x-shift of the gas-gap center w.r.t. quadruplet center
523 * Amg::getTranslateZ3D(m_offset) // z-shift to volume center (after m_delta!)
524 * Amg::getRotateY3D(-90 * CLHEP::deg) // x<->z because of GeoTrd definition
525 * Amg::getRotateZ3D(-90 * CLHEP::deg)); // x<->y for wires
526
527 m_surfaceData->m_layerCenters.emplace_back(m_surfaceData->m_layerTransforms.back().translation());
528 m_surfaceData->m_layerNormals.emplace_back(m_surfaceData->m_layerTransforms.back().linear() * Amg::Vector3D(0., 0., -1.));
529
530 //-------------------
531 // Strips
532 //-------------------
533
534 const double shift{layer%2 == 0 ? 0.01 : -0.01}; // 1st layer gets +0.01; layer numbering starts from 0 here!
535
536 // identifier of the first channel - strip plane
537 id = m_idHelper.channelID(getStationName(), getStationEta(), getStationPhi(), m_ml, layer + 1, 1, 1);
538
539 m_surfaceData->m_layerSurfaces.push_back(std::make_unique<Trk::PlaneSurface>(*this, id));
540
541 m_surfaceData->m_layerTransforms.push_back(absTransform() * m_delta * m_Xlg[layer] *Amg::Translation3D(shift, 0., m_offset)
542 *Amg::getRotateY3D(-90 * CLHEP::deg)); // x<->z because of GeoTrd definition
543
544 m_surfaceData->m_layerCenters.emplace_back(m_surfaceData->m_layerTransforms.back().translation());
545 m_surfaceData->m_layerNormals.emplace_back(m_surfaceData->m_layerTransforms.back().linear() * Amg::Vector3D(0., 0., -1.));
546
547 //-------------------
548 // Trigger Pads
549 //-------------------
550
551 // identifier of the first channel - pad plane
552 id = m_idHelper.channelID(getStationName(), getStationEta(), getStationPhi(), m_ml, layer + 1, 0, 1);
553
554 m_surfaceData->m_layerSurfaces.push_back(std::make_unique<Trk::PlaneSurface>(*this, id));
555
556 m_surfaceData->m_layerTransforms.push_back(absTransform() * m_delta * m_Xlg[layer] *
557 Amg::getTranslate3D(-shift, 0., m_offset)
558 * Amg::getRotateY3D(-90 * CLHEP::deg) // x<->z because of GeoTrd definition
559 * Amg::getRotateZ3D(-90 * CLHEP::deg)); // x<->y for pads
560
561 m_surfaceData->m_layerCenters.emplace_back(m_surfaceData->m_layerTransforms.back().translation());
562 m_surfaceData->m_layerNormals.emplace_back(m_surfaceData->m_layerTransforms.back().linear() * Amg::Vector3D(0., 0., -1.));
563 }
564 }
565
566
567 //============================================================================
569 if (m_idHelper.stationEta(id) != getStationEta()) return false;
570 if (m_idHelper.stationPhi(id) != getStationPhi()) return false;
571
572 if (m_idHelper.multilayerID(id) != m_ml) return false;
573
574 int gasgap = m_idHelper.gasGap(id);
575 if (gasgap < 1 || gasgap > m_nlayers) return false;
576
577 int strip = m_idHelper.channel(id);
578 if (strip < 1) return false;
579 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Strip && strip > m_etaDesign[gasgap - 1].nch) return false;
580 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Wire && strip > m_phiDesign[gasgap -1].nGroups) return false;
581 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Pad) {
582 const auto [etaId, phiId] = m_padDesign[gasgap -1].etaPhiId(strip);
583 if (etaId < 0 || phiId < 0) return false;
584 }
585 return true;
586 }
587
588
589 //============================================================================
591 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Pad) {
592 const MuonPadDesign* design = getPadDesign(id);
593 if (!design) {
594 ATH_MSG_WARNING( "no pad Design" );
595 return -1;
596 }
597 return design->channelWidth(Amg::Vector2D::Zero(), 0);
598 }
599
600 const MuonChannelDesign* design = getDesign(id);
601 if (!design) return -1;
602 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Strip) // sTGC strips
603 return design->inputPitch;
604 else if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Wire) // sTGC wires
605 return design->inputPitch * design->groupWidth; // wire Pitch * number of wires in a group
606 else
607 return -1;
608 }
609
610
611 //============================================================================
612 int sTgcReadoutElement::padNumber(const Amg::Vector2D& pos, const Identifier& id) const {
613 const MuonPadDesign* design = getPadDesign(id);
614 if (!design) {
615 ATH_MSG_WARNING( "no pad Design" );
616 return -1;
617 }
618 std::pair<int, int> pad(design->channelNumber(pos));
619 const sTgcIdHelper& id_helper{*manager()->stgcIdHelper()};
620 if (pad.first > 0 && pad.second > 0) {
621#ifndef NDEBUG
622 bool is_valid {true};
623#endif
624 const Identifier padID = id_helper.padID(id, id_helper.multilayer(id),
625 id_helper.gasGap(id), sTgcIdHelper::Pad, pad.first, pad.second
626#ifndef NDEBUG
627 , is_valid
628#endif
629
630 );
631 int channel = id_helper.channel(padID);
632 int padEta = id_helper.padEta(padID);
633 int padPhi = id_helper.padPhi(padID);
634 if (
635#ifndef NDEBUG
636 !is_valid ||
637#endif
638 padEta != pad.first || padPhi != pad.second) {
639
640 ATH_MSG_WARNING( " bad pad indices: input " << pad.first << " " << pad.second << " from ID " << padEta << " "
641 << padPhi );
642 return -1;
643 }
644 return channel;
645 }
646
647 ATH_MSG_WARNING(__LINE__<< " bad channelNumber" <<pad.first<<" "<<pad.second );
648
649 return -1;
650 }
651
652
653 //============================================================================
654 int sTgcReadoutElement::wireNumber(const Amg::Vector2D& pos, const Identifier& id) const {
655 const MuonChannelDesign* design = getDesign(id);
656 if (!design) {
657 ATH_MSG_WARNING( "no wire design when trying to get the wire number" );
658 return -1;
659 }
660 return design->wireNumber(pos);
661 }
662
663
664 //============================================================================
665 double sTgcReadoutElement::wirePitch(int gas_gap) const {
666 if (m_phiDesign.empty()) {
667 ATH_MSG_WARNING( "no wire design when trying to get the wire pitch" );
668 return -1.0;
669 }
670 return (m_phiDesign[gas_gap - 1]).inputPitch;
671 }
672
673
674 //============================================================================
676 double pos_wire = -9999.9;
677 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Wire) {
678 const MuonChannelDesign* design = getDesign(id);
679 if (!design) {
680 ATH_MSG_WARNING( "no wire design when trying to get the 1st wire position" );
681 return pos_wire;
682 }
683 pos_wire = design->firstPos();
684 } else {
685 ATH_MSG_WARNING( "attempt to retrieve the 1st wire position with a wrong identifier" );
686 }
687 return pos_wire;
688 }
689
690
691 //============================================================================
693 int nWires = -1;
694 if (m_idHelper.channelType(id) == sTgcIdHelper::sTgcChannelTypes::Wire) {
695 const MuonChannelDesign* design = getDesign(id);
696 if (!design) {
697 ATH_MSG_WARNING( "no wire design when trying to get the total number of wires" );
698 return nWires;
699 }
700 nWires = design->nch;
701 } else {
702 ATH_MSG_WARNING( "attempt to retrieve the number of wires with a wrong identifier" );
703 }
704 return nWires;
705 }
706
707
708 //============================================================================
710 int gg = m_idHelper.gasGap(id);
711 int channelType = m_idHelper.channelType(id);
712
713 // The assigned coordinate along the layer normal is at the center of the gas gap;
714 // wires are considered at x=0, while:
715 // for layers 1, 3 strips (pads) are shifted by +10 (-10) microns
716 // for layers 2, 4 strips (pads) are shifted by -10 (+10) microns
717 double shift{0.};
718 if (channelType != 2) shift = ((gg % 2) ^ (channelType==0)) ? 0.01 : -0.01;
719 const Amg::Vector3D locPos_ML = m_Xlg[gg - 1] * Amg::getTranslate3D(shift, 0., m_offset) * locPos;
720
721 ATH_MSG_DEBUG( "position coordinates in the gas-gap r.f.: " << Amg::toString(locPos) );
722 ATH_MSG_DEBUG( "position coordinates in the multilayer r.f.: " << Amg::toString(locPos_ML) );
723 return absTransform() * m_delta * locPos_ML;
724 }
725
726
727 //============================================================================
729 // amdb frame (s, z, t) = chamber frame (y, z, x)
730 if (aline) {
731 static const Amg::Transform3D permute{GeoTrf::GeoRotation{90.*Gaudi::Units::deg,90.*Gaudi::Units::deg, 0.}};
732 // The origin of the rotation axes is at the center of the active area
733 // in the z (radial) direction. Account for this shift in the definition
734 // of m_delta so that it can be applied on chamber frame coordinates.
735 m_ALinePar = &aline;
736 m_delta = Amg::getTranslateZ3D(m_offset)* permute*aline.delta()*
737 permute.inverse()*Amg::getTranslateZ3D(-m_offset);
738 ATH_MSG_DEBUG(idHelperSvc()->toStringDetEl(identify())<<" setup new alignment: "<<GeoTrf::toString(m_delta));
739 refreshCache();
740 } else {
742 }
743 }
744
745 //============================================================================
747 if (has_ALines()) {
748 m_ALinePar = nullptr;
749 m_delta = Amg::Transform3D::Identity();
750 refreshCache();
751 }
752 }
753
754 //============================================================================
756 ATH_MSG_DEBUG("Setting B-line for " <<idHelperSvc()->toStringDetEl(identify())<<" "<<bLine);
757 m_BLinePar = &bLine;
758 }
759
760 //============================================================================
762
763 // note: amdb frame (s, z, t) = chamber frame (y, z, x)
764 if (!has_BLines()) return;
765
766 double t0 = locPosML.x();
767 double s0 = locPosML.y();
768 double z0 = locPosML.z();
769 double width = getSsize() + (getLongSsize() - getSsize())*(z0/getRsize() + 0.5); // because z0 is in [-length/2, length/2]
770
771 double s_rel = s0/(width/2.); // in [-1, 1]
772 double z_rel = z0/(getRsize()/2.); // in [-1, 1]
773 double t_rel = t0/(getZsize()/2.); // in [-1, 1]
774
775 // b-line parameters
776 using Parameter = BLinePar::Parameter;
777 const double bp = m_BLinePar->getParameter(Parameter::bp);
778 const double bn = m_BLinePar->getParameter(Parameter::bn);
779 const double sp = m_BLinePar->getParameter(Parameter::sp);
780 const double sn = m_BLinePar->getParameter(Parameter::sn);
781 const double tw = m_BLinePar->getParameter(Parameter::tw);
782 const double eg = m_BLinePar->getParameter(Parameter::eg)*1.e-3;
783 const double ep = m_BLinePar->getParameter(Parameter::ep)*1.e-3;
784 const double en = m_BLinePar->getParameter(Parameter::en)*1.e-3;
785
786 double ds{0.}, dz{0.}, dt{0.};
787
788 if (bp != 0 || bn != 0)
789 dt += 0.5*(s_rel*s_rel - 1)*((bp + bn) + (bp - bn)*z_rel);
790
791 if (sp != 0 || sn != 0)
792 dt += 0.5*(z_rel*z_rel - 1)*((sp + sn) + (sp - sn)*s_rel);
793
794 if (tw != 0) {
795 dt -= tw*s_rel*z_rel;
796 dz += tw*s_rel*t_rel*getZsize()/getRsize();
797 }
798
799 if (eg != 0) {
800 dt += t0*eg;
801 ds += s0*eg;
802 dz += z0*eg;
803 }
804
805 if (ep != 0 || en != 0) {
806 // the formulas below differ from those in Christoph's talk
807 // because are origin for NSW is at the center of the chamber,
808 // whereas in the talk (i.e. MDTs), it is at the bottom!
809 double delta = s_rel*s_rel * ((ep + en)*s_rel/6 + (ep - en)/4);
810 double phi = s_rel * ((ep + en)*s_rel + (ep - en)) / 2;
811 dt += phi*t0;
812 ds += delta*width/2;
813 dz += phi*z0;
814 }
815
816 locPosML[0] += dt;
817 locPosML[1] += ds;
818 locPosML[2] += dz;
819 }
820
821
822 //============================================================================
823 void sTgcReadoutElement::spacePointPosition(const Identifier& layerId, double locXpos, double locYpos, Amg::Vector3D& pos) const {
824
825 pos = Amg::Vector3D(locXpos, locYpos, 0.);
826
827 const MuonChannelDesign* design = getDesign(layerId);
828 if (!design) {
829 ATH_MSG_WARNING( "Unable to get MuonChannelDesign, therefore cannot provide position corrections. Returning." );
830 return;
831 }
832
833#ifndef SIMULATIONBASE
834 //*********************
835 // As-Built (MuonNswAsBuilt is not included in AthSimulation)
836 //*********************
837 if(manager()->getsTGCAsBuilt() && design->type == MuonChannelDesign::ChannelType::etaStrip){
838#if __GNUC__ >= 13
839// Avoid a warning seen with -march=x86-64-v3.
840// This has been cleaned up in eigen after 3.4.0.
841# pragma GCC diagnostic push
842# pragma GCC diagnostic ignored "-Warray-bounds"
843#endif
844 pos.head(2) = manager()->getsTGCAsBuilt()->correctPosition(layerId, pos.head(2));
845#if __GNUC__ >= 13
846# pragma GCC diagnostic pop
847#endif
848 }
849#ifndef NDEBUG
850 else {
851 MsgStream log(Athena::getMessageSvc(), "sTgcReadoutElement");
852 if (log.level() <= MSG::DEBUG) {
853 log << MSG::DEBUG << "No as-built corrections provided for stEta: "<<getStationEta() << " stPhi: "<<getStationPhi()<<" ml: "<<m_ml<< endmsg;
854 }
855 }
856#endif
857#endif
858
859
860 //*********************
861 // B-Lines
862 //*********************
863 if (has_BLines()) {
864 // go to the muultilayer frame
865 Amg::Transform3D trfToML = m_delta.inverse()*absTransform().inverse()*transform(layerId);
866 pos = trfToML*pos;
867 posOnDefChamber(pos);
868 // back to the layer reference frame from where we started
869 pos = trfToML.inverse()*pos;
870
871 }
872 }
873
874} // namespace MuonGM
Scalar phi() const
phi method
#define endmsg
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
Definition of the abstract IRDBAccessSvc interface.
std::shared_ptr< IRDBRecordset > IRDBRecordset_ptr
Definition of the abstract IRDBRecord interface.
Definition of the abstract IRDBRecordset interface.
static Double_t sp
static Double_t s0
static Double_t t0
const double width
const std::string & GetName() const
Amg::Transform3D delta() const
Returns the final transformations of the A lines.
Definition ALinePar.cxx:36
IRDBRecord is one record in the IRDBRecordset object.
Definition IRDBRecord.h:28
virtual double getDouble(std::string_view fieldName) const =0
Get double field value.
virtual const std::string & getString(std::string_view fieldName) const =0
Get string field value.
virtual unsigned int size() const =0
std::unique_ptr< SurfaceData > m_surfaceData
MuonClusterReadoutElement(GeoVFullPhysVol *pv, MuonDetectorManager *mgr, Trk::DetectorElemType detType)
virtual const Amg::Transform3D & transform() const override
Return local to global transform.
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
void setIdentifier(const Identifier &id)
Sets the Identifier, hashes & station names.
Identifier identify() const override final
Returns the ATLAS Identifier of the MuonReadOutElement.
const MuonPadDesign * getPadDesign(const Identifier &id) const
returns the MuonChannelDesign class for the given identifier
virtual void fillCache() override final
function to fill tracking cache
const MuonChannelDesign * getDesign(const Identifier &id) const
returns the MuonChannelDesign class for the given identifier
int numberOfWires(const Identifier &id) const
Get the total number of wires (single wires) of a chamber.
double channelPitch(const Identifier &id) const
Channel pitch.
sTgcReadoutElement(GeoVFullPhysVol *pv, const std::string &stName, int zi, int fi, int mL, MuonDetectorManager *mgr)
constructor
Amg::Vector3D localToGlobalCoords(const Amg::Vector3D &locPos, Identifier id) const
simHit local (SD) To Global position - to be used by MuonGeoAdaprors only
void setBLinePar(const BLinePar &bLine)
read B-line (chamber-deformation) parameters
void initDesign(double thickness)
initialize the design classes for this readout element
virtual bool spacePointPosition(const Identifier &phiId, const Identifier &etaId, Amg::Vector2D &pos) const override final
space point position for a given pair of phi and eta identifiers The LocalPosition is expressed in th...
double positionFirstWire(const Identifier &id) const
Get the local position of the first wire of the chamber corresponding to the identifier.
virtual bool containsId(const Identifier &id) const override final
function to be used to check whether a given Identifier is contained in the readout element
void posOnDefChamber(Amg::Vector3D &locPosML) const
transform a position (in chamber-frame coordinates) to the deformed-chamber geometry
int wireNumber(const Amg::Vector2D &pos, const Identifier &id) const
wire number corresponding to local position
void setChamberLayer(int ml)
set methods only to be used by MuonGeoModel
int padNumber(const Amg::Vector2D &pos, const Identifier &id) const
pad number corresponding to local position
Definition node.h:24
Amg::Vector2D correctPosition(const Identifier &channelId, const Amg::Vector2D &pos) const
sTGCReadoutParameters & GetReadoutParameters()
MuonGM::sTGC_Technology * GetTechnology()
sTGCDetectorDescription * Get_sTGCDetector(char type, int ieta, int iphi, int layer=1, char side='A')
int padPhi(const Identifier &id) const
int multilayer(const Identifier &id) const
int padEta(const Identifier &id) const
int channel(const Identifier &id) const override
Identifier padID(int stationName, int stationEta, int stationPhi, int multilayer, int gasGap, int channelType, int padEta, int padPhi) const
int gasGap(const Identifier &id) const override
get the hashes
Amg::Transform3D getTranslate3D(const double X, const double Y, const double Z)
: Returns a shift transformation along an arbitrary axis
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Amg::Transform3D getTranslateZ3D(const double Z)
: Returns a shift transformation along the z-axis
Amg::Transform3D getRotateZ3D(double angle)
get a rotation transformation around Z-axis
Eigen::Affine3d Transform3D
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
Amg::Transform3D getRotateY3D(double angle)
get a rotation transformation around Y-axis
Eigen::Translation< double, 3 > Translation3D
IMessageSvc * getMessageSvc(bool quiet=false)
Ensure that the Athena extensions are properly loaded.
Definition GeoMuonHits.h:27
int nStrips(const MuonGM::TgcReadoutElement &readoutEle, int layer)
Ensure that the ATLAS eigen extensions are properly loaded.
Definition index.py:1
int wireNumber(const Amg::Vector2D &pos) const
calculate the sTGC wire number. The method can return a value outside the range [1,...
double firstPos() const
Returns the position of the first strip along the x-axis.
Parameters defining the design of the readout sTGC pads.
std::pair< int, int > channelNumber(const Amg::Vector2D &pos) const
calculate local channel number, range 1=nstrips like identifiers.
double channelWidth(const Amg::Vector2D &pos, bool measPhi, bool preciseMeas=false) const
calculate local channel width
std::vector< int > firstWireGroup
std::vector< double > firstPadPhiDivision_A
std::vector< double > firstStripWidth
std::vector< double > firstPadH
std::vector< double > firstWire
std::vector< double > padH
std::vector< double > wireCutout
std::vector< double > PadPhiShift_A
std::vector< double > nPadH