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