ATLAS Offline Software
Loading...
Searching...
No Matches
TgcDigitMaker.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#include "TgcDigitMaker.h"
6
7#include <cmath>
8#include <fstream>
9#include <iostream>
10#include <memory>
11
12#include "AthenaBaseComps/AthCheckMacros.h"
13#include "CLHEP/Random/RandFlat.h"
14#include "CLHEP/Random/RandomEngine.h"
15#include "CLHEP/Units/SystemOfUnits.h"
16#include "GaudiKernel/MsgStream.h"
17#include "MuonDigitContainer/TgcDigitCollection.h"
18#include "MuonIdHelpers/TgcIdHelper.h"
19#include "MuonReadoutGeometry/MuonDetectorManager.h"
20#include "MuonReadoutGeometry/TgcReadoutElement.h"
21#include "MuonSimEvent/TGCSimHit.h"
22#include "MuonSimEvent/TgcHitIdHelper.h"
23#include "PathResolver/PathResolver.h"
24
25// run number from geometry DB
26#include "GaudiKernel/Bootstrap.h"
27#include "GeoModelInterfaces/IGeoModelSvc.h"
28#include "RDBAccessSvc/IRDBAccessSvc.h"
29#include "RDBAccessSvc/IRDBRecord.h"
30#include "RDBAccessSvc/IRDBRecordset.h"
31
32//---------------------------------------------------
33// Constructor and Destructor
34//---------------------------------------------------
35
36//----- Constructor
37TgcDigitMaker::TgcDigitMaker(const TgcHitIdHelper* hitIdHelper,
38 const MuonGM::MuonDetectorManager* mdManager, const bool doFourBunch)
39 : AthMessaging("TgcDigitMaker"), m_doFourBunchDigitization(doFourBunch) {
40 m_hitIdHelper = hitIdHelper;
41 m_mdManager = mdManager;
42 m_idHelper = nullptr;
43 m_efficiency[kWIRE] = m_efficiency[kSTRIP] =
44 1.000; // 100% efficiency for TGCSimHit_p1
45 m_gateTimeWindow[kOUTER][kWIRE] =
46 29.32; // 29.32ns = 26ns + 4 * 0.83ns(outer station)
47 m_gateTimeWindow[kOUTER][kSTRIP] =
48 40.94; // 40.94ns = 26ns + 18 * 0.83ns(outer station)
49 m_gateTimeWindow[kINNER][kWIRE] =
50 33.47; // 33.47ns = 26ns + 9 * 0.83ns(inner station)
51 m_gateTimeWindow[kINNER][kSTRIP] =
52 45.09; // 45.09ns = 26ns + 23 * 0.83ns(inner station)
53 m_bunchCrossingTime = 24.95; // 24.95 ns =(40.08 MHz)^(-1)
54}
55
56//------------------------------------------------------
57// Initialize
58//------------------------------------------------------
59StatusCode TgcDigitMaker::initialize() {
60 // Initialize TgcIdHelper
61 if (!m_hitIdHelper) {
62 m_hitIdHelper = TgcHitIdHelper::GetHelper();
63 }
64
65 // initialize the TGC identifier helper
66 m_idHelper = m_mdManager->tgcIdHelper();
67
68 ATH_CHECK(readFileOfTimeJitter());
69
70 // Read share/TGC_Digitization_energyThreshold.dat file and store values in
71 // m_energyThreshold.
72 ATH_CHECK(readFileOfEnergyThreshold());
73
74 // Read share/TGC_Digitization_deadChamber.dat file and store values in
75 // m_isDeadChamber.
76 ATH_CHECK(readFileOfDeadChamber());
77
78 // Read share/TGC_Digitization_StripPosition.dat file and store values in
79 // m_StripPosition.
80 ATH_CHECK(readFileOfStripPosition());
81
82 return StatusCode::SUCCESS;
83}
84
85//---------------------------------------------------
86// Execute Digitization
87//---------------------------------------------------
88TgcDigitCollection* TgcDigitMaker::executeDigi(
89 const TGCSimHit* hit, const double globalHitTime,
90 const TgcDigitASDposData* ASDpos, const TgcDigitTimeOffsetData* TOffset,
91 const TgcDigitCrosstalkData* Crosstalk,
92 CLHEP::HepRandomEngine* rndmEngine) {
93 // timing constant parameters
94 constexpr float sensor_propagation_time =
95 3.3 * CLHEP::ns /
96 CLHEP::m; // Until MC10, 8.5*ns/m for wire, 8.7*ns/m for strip.
97 // Since MC11, 3.3*ns/m (the speed of light) is used
98 // from the Z->mumu data/MC comparison.
99 constexpr float cable_propagation_time = 5.0 * CLHEP::ns / CLHEP::m;
100
101 // position constant parameters
102 constexpr float wire_pitch = 1.8 * CLHEP::mm;
103 constexpr float zwidth_frame = 17. * CLHEP::mm;
104
105 ATH_MSG_DEBUG("executeDigi() Got HIT Id.");
106
108 int Id = hit->TGCid();
109 std::string stationName = m_hitIdHelper->GetStationName(Id);
110 int stationEta = m_hitIdHelper->GetStationEta(Id);
111 int stationPhi = m_hitIdHelper->GetStationPhi(Id);
112 int ilyr = m_hitIdHelper->GetGasGap(Id);
113
114 // Check the chamber is dead or not.
115 if (isDeadChamber(stationName, stationEta, stationPhi, ilyr))
116 return nullptr;
117
118 const Identifier elemId =
119 m_idHelper->elementID(stationName, stationEta, stationPhi);
120 ATH_MSG_DEBUG("executeDigi() - element identifier is: "
121 << m_idHelper->show_to_string(elemId));
122
123 const MuonGM::TgcReadoutElement* tgcChamber =
124 m_mdManager->getTgcReadoutElement(elemId);
125 if (!tgcChamber) {
126 ATH_MSG_WARNING("executeDigi() - no ReadoutElement found for "
127 << m_idHelper->show_to_string(elemId));
128 return nullptr;
129 }
130
131 IdContext tgcContext = m_idHelper->module_context();
132 IdentifierHash coll_hash;
133 if (m_idHelper->get_hash(elemId, coll_hash, &tgcContext)) {
134 ATH_MSG_WARNING("Unable to get TGC hash id from TGC Digit collection "
135 << "context begin_index = " << tgcContext.begin_index()
136 << " context end_index = " << tgcContext.end_index()
137 << " the identifier is " << elemId);
138 }
139
140 std::unique_ptr<TgcDigitCollection> digits =
141 std::make_unique<TgcDigitCollection>(elemId, coll_hash);
142
143 const Amg::Vector3D centreChamber = tgcChamber->globalPosition();
144 float height = tgcChamber->getRsize();
145 float hmin = sqrt(pow(centreChamber.x(), 2) + pow(centreChamber.y(), 2)) -
146 height / 2.;
147
148 // Direction cosine of incident angle of this track
149 Amg::Vector3D direCos = hit->localDireCos();
150
151 // local position
152 Amg::Vector3D localPos = hit->localPosition();
153
154 // Local z direction is global r direction.
155 float distanceZ = 1.4 * CLHEP::mm / direCos[0] * direCos[2];
156 float zLocal = localPos.z() + distanceZ;
157
158 // Local y direction is global phi direction.
159 float distanceY = 1.4 * CLHEP::mm / direCos[0] * direCos[1];
160 // This ySign depends on the implementation of TGC geometry in G4 simulation
161 // left-handed coordinate in A side(+z, stationEta>0)
162 float ySign = (stationEta < 0) ? +1. : -1.;
163 float yLocal = ySign * (localPos.y() + distanceY);
164
165 // Time of flight correction for each chamber
166 // the offset is set to 0 for ASD located at the postion where tof is
167 // minimum in each chamber, i.e. the ASD at the smallest R in each chamber
168 double tofCorrection =
169 (sqrt(pow(hmin + zwidth_frame, 2) + pow(centreChamber.z(), 2)) /
170 (299792458. * (CLHEP::m / CLHEP::s))); // FIXME use CLHEP::c_light
171
172 // bunch time
173 float bunchTime = globalHitTime - tofCorrection;
174
175 static const float jitterInitial = 9999.;
176 float jitter = jitterInitial; // calculated at wire group calculation and
177 // used also for strip calculation
178
179 int iStationName = getIStationName(stationName);
180
181 double energyDeposit = hit->energyDeposit(); // Energy deposit in MeV
182 // If TGCSimHit_p1 is used, energyDeposit is the default value of -9999.
183 // If TGCSimHit_p2 is used, energyDeposit is equal to or greater than 0.
184 // Therefore, the energyDeposit value can be used to distinguish
185 // TGCSimHit_p1 and TGCSimHit_p2. For TGCSimHit_p1, old efficiency check
186 // with only isStrip variable is used. For TGCSimHit_p2, new efficiency
187 // check with chamber dependent energy threshold is used.
188
190
191 if ((energyDeposit >= -1. &&
192 efficiencyCheck(
193 stationName, stationEta, stationPhi, ilyr, kWIRE,
194 energyDeposit)) || // efficiency check for TGCSimHit_p2 at first,
195 (energyDeposit < -1. &&
196 efficiencyCheck(
197 kWIRE, rndmEngine))) { // Old efficiencyCheck for TGCSimHit_p1
198
199 int iWireGroup[2];
200 float posInWireGroup[2] = {0., 0.};
201 for (int iPosition = 0; iPosition < 2; iPosition++) {
202 int nWireOffset = (std::abs(stationEta) == 5 ||
203 stationName.compare(0, 2, "T4") == 0)
204 ? 1
205 : 0;
206 // for chambers in which only the first wire is not connected : 1
207 // for chambers in which the first and last wires are not connected
208 // OR all wires are connected : 0
209
210 double zPosInSensArea =
211 zLocal + static_cast<double>(tgcChamber->nWires(ilyr) -
212 nWireOffset) *
213 wire_pitch / 2.;
214
215 // check a hit in the sensitive area
216 if (zPosInSensArea < 0. ||
217 zPosInSensArea > tgcChamber->nWires(ilyr) * wire_pitch) {
218 iWireGroup[iPosition] = 0;
219 posInWireGroup[iPosition] = 0.;
220 ATH_MSG_DEBUG(
221 "executeDigi(): Wire: Hit position located at outside of a "
222 "sensitive volume, "
223 << " id: " << stationName << "/" << stationEta << "/"
224 << stationPhi << "/" << ilyr
225 << " position: " << zPosInSensArea << " xlocal: " << zLocal
226 << " dir cosine: " << direCos[0] << "/" << direCos[1] << "/"
227 << direCos[2]
228 << " active region: " << height - zwidth_frame * 2.);
229 } else {
230 int igang = 1;
231 int wire_index = 0;
232 while (wire_pitch * (static_cast<float>(wire_index)) <
233 zPosInSensArea &&
234 igang <= tgcChamber->nWireGangs(ilyr)) {
235 wire_index += tgcChamber->nWires(ilyr, igang);
236 igang++;
237 }
238 posInWireGroup[iPosition] =
239 (zPosInSensArea / wire_pitch -
240 (static_cast<float>(wire_index))) /
241 (static_cast<float>(
242 tgcChamber->nWires(ilyr, igang - 1))) +
243 1.;
244
245 iWireGroup[iPosition] = ((1 == igang) ? 1 : igang - 1);
246 }
247 }
248
249 unsigned int jWG[2] = {
250 (iWireGroup[0] <= iWireGroup[1]) ? (unsigned int)(0)
251 : (unsigned int)(1),
252 (iWireGroup[0] <= iWireGroup[1]) ? (unsigned int)(1)
253 : (unsigned int)(0)};
254 int iWG[2] = {iWireGroup[jWG[0]], iWireGroup[jWG[1]]};
255 float posInWG[2] = {posInWireGroup[jWG[0]], posInWireGroup[jWG[1]]};
256 if (iWG[0] != iWG[1]) {
257 ATH_MSG_DEBUG("executeDigi(): Multihits found in WIRE GROUPs:"
258 << iWG[0] << " " << iWG[1]);
259 }
260
261 // === BC tagging from the hit timing ===
262 for (int iwg = iWG[0]; iwg <= iWG[1]; iwg++) {
263 if (1 <= iwg && iwg <= tgcChamber->nWireGangs(ilyr)) {
264 // timing window offset
265 float wire_timeOffset =
266 (TOffset != nullptr)
267 ? this->getTimeOffset(TOffset, iStationName, stationEta,
268 kWIRE)
269 : 0.;
270 // EI/FI has different offset between two ASDs.
271 if (iStationName > 46)
272 wire_timeOffset +=
273 (iwg < 17) ? 0.5 * CLHEP::ns : -0.5 * CLHEP::ns;
274
275 // TGC response time calculation
276 float jit = timeJitter(direCos, rndmEngine);
277 if (jit < jitter)
278 jitter = jit;
279 float ySignPhi = (stationPhi % 2 == 1) ? -1. : +1.;
280 // stationPhi%2==0 : +1. : ASD attached at the smaller phi side
281 // of TGC stationPhi%2==1 : -1. : ASD attached at the larger phi
282 // side of TGC
283 float wTimeDiffByRadiusOfInner =
284 this->timeDiffByCableRadiusOfInner(iStationName, stationPhi,
285 iwg);
286 double digit_time =
287 bunchTime + jitter + wTimeDiffByRadiusOfInner;
288 float wASDDis{-1000.};
289 if (ASDpos != nullptr) {
290 wASDDis = this->getDistanceToAsdFromSensor(
291 ASDpos, iStationName, abs(stationEta), stationPhi,
292 kWIRE, iwg, zLocal);
293 float wCableDis =
294 wASDDis + (ySignPhi * yLocal +
295 tgcChamber->chamberWidth(zLocal) / 2.);
296 float wPropTime =
297 sensor_propagation_time *
298 (ySignPhi * yLocal +
299 tgcChamber->chamberWidth(zLocal) / 2.) +
300 cable_propagation_time * wASDDis;
301 digit_time +=
302 this->getSigPropTimeDelay(wCableDis) + wPropTime;
303 }
304
305 TgcStation station =
306 (m_idHelper->stationName(elemId) > 46) ? kINNER : kOUTER;
307 uint16_t bctag =
308 bcTagging(digit_time, m_gateTimeWindow[station][kWIRE],
309 wire_timeOffset);
310
311 if (bctag == 0x0) {
312 ATH_MSG_DEBUG(
313 "WireGroup: digitized time "
314 << digit_time << " ns is outside of time window from "
315 << wire_timeOffset << ". bunchTime: " << bunchTime
316 << " time jitter: " << jitter << " propagation time: "
317 << sensor_propagation_time *
318 (ySignPhi * yLocal +
319 tgcChamber->chamberWidth(zLocal) / 2.) +
320 cable_propagation_time * wASDDis
321 << " time difference by cable radius of inner station: "
322 << wTimeDiffByRadiusOfInner);
323 } else {
324 Identifier newId = m_idHelper->channelID(
325 stationName, stationEta, stationPhi, ilyr, (int)kWIRE,
326 iwg);
327 addDigit(newId, bctag, digits.get());
328
329 if (iwg == iWG[0]) {
330 randomCrossTalk(Crosstalk, elemId, ilyr, kWIRE, iwg,
331 posInWG[0], digit_time, wire_timeOffset,
332 rndmEngine, digits.get());
333 }
334
335 ATH_MSG_DEBUG(
336 "WireGroup: newid breakdown digitTime x/y/z direcos "
337 "height_gang bctag: "
338 << newId << " " << stationName << "/" << stationEta
339 << "/" << stationPhi << "/" << ilyr << "/"
340 << "WIRE/" << iwg << " " << digit_time << " "
341 << localPos.x() << "/" << localPos.y() << "/"
342 << localPos.z() << " " << direCos.x() << "/"
343 << direCos.y() << "/" << direCos.z() << " " << height
344 << " " << tgcChamber->nWires(ilyr, iwg) << " "
345 << bctag);
346 }
347 } else {
348 ATH_MSG_DEBUG(
349 "Wire Gang id is out of range. id, x/y/z, height: "
350 << iwg << " " << localPos.x() << "/" << localPos.y() << "/"
351 << localPos.z() << " " << height);
352 }
353 }
354 } // end of wire group calculation
355
357 TgcSensor sensor = kSTRIP;
358
359 if ((ilyr != 2 || (stationName.compare(0, 2, "T1") !=
360 0)) && // no stip in middle layers of T1*
361 ((energyDeposit < -1. &&
362 efficiencyCheck(
363 sensor, rndmEngine)) || // Old efficiencyCheck for TGCSimHit_p1.
364 (energyDeposit >= -1. &&
365 efficiencyCheck(
366 stationName, stationEta, stationPhi, ilyr, sensor,
367 energyDeposit))) // New efficiencyCheck for TGCSimHit_p2
368 ) {
369
370
371 int iStrip[2];
372 float posInStrip[2] = {0., 0.};
373 // Take into account of charge spread on cathod plane
374 for (int iPosition = 0; iPosition < 2; iPosition++) {
375
376 // check a hit in the sensitive area
377 float zPos = zLocal + height / 2. - zwidth_frame;
378 if (zPos < 0.) {
379 iStrip[iPosition] = 0;
380 posInStrip[iPosition] = 0.;
381 ATH_MSG_DEBUG(
382 "Strip: Hit position located at outside of a sensitive "
383 "volume, id, position, xlocal0/1, dir cosine: "
384 << stationName << "/" << stationEta << "/" << stationPhi
385 << "/" << ilyr << zPos << " " << zLocal << " " << direCos[0]
386 << "/" << direCos[1] << "/" << direCos[2]);
387 } else if (zPos > height - zwidth_frame * 2.) {
388 iStrip[iPosition] = tgcChamber->nStrips(ilyr) + 1;
389 posInStrip[iPosition] = 0.;
390 ATH_MSG_DEBUG(
391 "Strip: Hit position located at outside of a sensitive "
392 "volume, id, position, active region: "
393 << stationName << "/" << stationEta << "/" << stationPhi
394 << "/" << ilyr << zPos << " "
395 << height - zwidth_frame * 2.);
396 } else { // sensitive area
397
398 //
399 // for layout P03
400 //
401 // number of strips in exclusive phi coverage of a chamber in
402 // T[1-3] and T4
403 const float dphi = tgcChamber->stripDeltaPhi();
404 float phiLocal = atan2(yLocal, zLocal + height / 2. + hmin);
405
406 ATH_MSG_DEBUG(
407 "dphi, phiLocal, yLocal, zLocall+ height/2.+hmin: "
408 << dphi << " " << phiLocal << " " << yLocal << " "
409 << zLocal + height / 2. + hmin);
410
411 int istr = 0;
412 if ((stationEta > 0 && ilyr == 1) ||
413 (stationEta < 0 && ilyr != 1)) {
414 istr = static_cast<int>(floor(phiLocal / dphi + 15.75)) + 1;
415 posInStrip[iPosition] =
416 phiLocal / dphi + 15.75 - static_cast<float>(istr - 1);
417 if (istr > 30) { // treatment for two half strips
418 istr = static_cast<int>(floor(
419 (phiLocal - dphi * 14.25) / (dphi / 2.))) +
420 31;
421 posInStrip[iPosition] =
422 (phiLocal - dphi * 14.25) / (dphi / 2.) -
423 static_cast<float>(istr - 31);
424 }
425 } else {
426 istr = static_cast<int>(floor(phiLocal / dphi + 16.25)) + 1;
427 posInStrip[iPosition] =
428 phiLocal / dphi + 16.25 - static_cast<float>(istr - 1);
429 if (istr < 3) { // treatment for two half strips
430 istr = static_cast<int>(floor(
431 (phiLocal + dphi * 14.25) / (dphi / 2.))) +
432 3;
433 posInStrip[iPosition] =
434 (phiLocal + dphi * 14.25) / (dphi / 2.) -
435 static_cast<float>(istr - 3);
436 }
437 }
438 if (istr < 1) {
439 istr = 0;
440 posInStrip[iPosition] = 0.;
441 } else if (32 < istr) {
442 istr = 33;
443 posInStrip[iPosition] = 0.;
444 }
445 iStrip[iPosition] = istr;
446 } // sensitive area
447 }
448
449 unsigned int jStr[2] = {
450 (iStrip[0] <= iStrip[1]) ? (unsigned int)(0) : (unsigned int)(1),
451 (iStrip[0] <= iStrip[1]) ? (unsigned int)(1) : (unsigned int)(0)};
452 int iStr[2] = {iStrip[jStr[0]], iStrip[jStr[1]]};
453 float posInStr[2] = {posInStrip[jStr[0]], posInStrip[jStr[1]]};
454 if (iStr[0] != iStr[1]) {
455 ATH_MSG_DEBUG("Multihits found in STRIPs:" << iStr[0] << " "
456 << iStr[1]);
457 }
458
459 // BC tagging from the timing
460 float strip_timeOffset =
461 (TOffset != nullptr)
462 ? this->getTimeOffset(TOffset, iStationName, stationEta, kSTRIP)
463 : 0.;
464
465 for (int istr = iStr[0]; istr <= iStr[1]; istr++) {
466 if (1 <= istr && istr <= 32) {
467 // TGC response time calculation
468 if (jitter > jitterInitial - 0.1) {
469 jitter = timeJitter(direCos, rndmEngine);
470 }
471 float zSignEta =
472 (abs(stationEta) % 2 == 1 && abs(stationEta) != 5) ? -1.
473 : +1.;
474 // if(abs(stationEta)%2 == 1 && abs(stationEta) != 5) : -1. :
475 // ASD attached at the longer base of TGC else : +1. : ASD
476 // attached at the shorter base of TGC
477 float sTimeDiffByRadiusOfInner =
478 this->timeDiffByCableRadiusOfInner(iStationName, stationPhi,
479 istr);
480 float sDigitTime =
481 bunchTime + jitter +
482 sensor_propagation_time *
483 (height / 2. + zwidth_frame + zSignEta * zLocal) +
484 sTimeDiffByRadiusOfInner;
485 float sASDDis{-1000};
486 if (ASDpos != nullptr) {
487 sASDDis = this->getDistanceToAsdFromSensor(
488 ASDpos, iStationName, abs(stationEta), stationPhi,
489 sensor, istr,
490 getStripPosition(stationName, abs(stationEta), istr));
491 float sCableDis = sASDDis + (height / 2. + zwidth_frame +
492 zSignEta * zLocal);
493 sDigitTime += this->getSigPropTimeDelay(sCableDis) +
494 sASDDis * cable_propagation_time;
495 }
496
497 TgcStation station =
498 (m_idHelper->stationName(elemId) > 46) ? kINNER : kOUTER;
499 uint16_t bctag =
500 bcTagging(sDigitTime, m_gateTimeWindow[station][sensor],
501 strip_timeOffset);
502
503 if (bctag == 0x0) {
504 ATH_MSG_DEBUG(
505 "Strip: Digitized time is outside of time window. "
506 << sDigitTime << " bunchTime: " << bunchTime
507 << " time jitter: " << jitter << " propagation time: "
508 << sensor_propagation_time *
509 (height / 2. + zwidth_frame + zSignEta * zLocal)
510 << " time difference by cable radius of inner station: "
511 << sTimeDiffByRadiusOfInner);
512 } else {
513 Identifier newId = m_idHelper->channelID(
514 stationName, stationEta, stationPhi, ilyr, (int)sensor,
515 istr);
516 addDigit(newId, bctag, digits.get());
517
518 if (istr == iStr[0]) {
519 randomCrossTalk(Crosstalk, elemId, ilyr, sensor,
520 iStr[0], posInStr[0], sDigitTime,
521 strip_timeOffset, rndmEngine,
522 digits.get());
523 }
524
525 ATH_MSG_DEBUG(
526 "Strip: newid breakdown digitTime x/y/z direcos "
527 "r_center bctag: "
528 << newId << " " << stationName << "/" << stationEta
529 << "/" << stationPhi << "/" << ilyr << "/" << sensor
530 << "/" << istr << " " << sDigitTime << " "
531 << localPos.x() << "/" << localPos.y() << "/"
532 << localPos.z() << " " << direCos.x() << "/"
533 << direCos.y() << "/" << direCos.z() << " "
534 << height / 2. + hmin << " " << bctag);
535 }
536 } else {
537 ATH_MSG_DEBUG("Strip id is out of range: gap, id: "
538 << ilyr << " " << istr);
539 }
540 }
541 } // end of strip number calculation
542
543 return digits.release();
544}
545
546//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
547StatusCode TgcDigitMaker::readFileOfTimeJitter() {
548 const char* const fileName = "TGC_Digitization_timejitter.dat";
549 std::string fileWithPath = PathResolver::find_file(fileName, "DATAPATH");
550
551 std::ifstream ifs;
552 if (!fileWithPath.empty()) {
553 ifs.open(fileWithPath.c_str(), std::ios::in);
554 } else {
555 ATH_MSG_FATAL("readFileOfTimeJitter(): Could not find file "
556 << fileName);
557 return StatusCode::FAILURE;
558 }
559
560 if (ifs.bad()) {
561 ATH_MSG_FATAL("readFileOfTimeJitter(): Could not open file "
562 << fileName);
563 return StatusCode::FAILURE;
564 }
565
566 int angle = 0;
567 int bins = 0;
568 int i = 0;
569 float prob = 0.;
570 bool verbose = msgLvl(MSG::VERBOSE);
571
572 while (ifs.good()) {
573 ifs >> angle >> bins;
574 if (ifs.eof())
575 break;
576 ATH_MSG_VERBOSE(
577 "readFileOfTimeJitter(): Timejitter, angle, Number of bins, prob. "
578 "dist.: "
579 << angle << " " << bins << " ");
580 m_vecAngle_Time.resize(i + 1);
581 for (int j = 0; j < 41 /*bins*/; j++) {
582 ifs >> prob;
583 m_vecAngle_Time[i].push_back(prob);
584 if (j == 0 && verbose)
585 msg(MSG::VERBOSE) << "readFileOfTimeJitter(): ";
586 if (verbose)
587 msg(MSG::VERBOSE) << prob << " ";
588 }
589 if (verbose)
590 msg(MSG::VERBOSE) << endmsg;
591 i++;
592 }
593 ifs.close();
594 return StatusCode::SUCCESS;
595}
596//+++++++++++++++++++++++++++++++++++++++++++++++
597float TgcDigitMaker::timeJitter(const Amg::Vector3D& direCosLocal,
598 CLHEP::HepRandomEngine* rndmEngine) const {
599 float injectionAngle =
600 atan2(fabs(direCosLocal[2]), fabs(direCosLocal[0])) / CLHEP::degree;
601
602 int ithAngle = static_cast<int>(injectionAngle / 5.);
603 float wAngle = injectionAngle / 5. - static_cast<float>(ithAngle);
604 int jthAngle;
605 if (ithAngle > 11) {
606 ithAngle = 12;
607 jthAngle = 12;
608 } else {
609 jthAngle = ithAngle + 1;
610 }
611
612 float jitter = 0.;
613 float prob = 1.;
614 float probRef = 0.;
615
616 while (prob > probRef) {
617 prob = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0);
618 jitter = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) * 40. *
619 CLHEP::ns; // trial time jitter
620 int ithJitter = static_cast<int>(jitter);
621 // probability distribution calculated from weighted sum between
622 // neighboring bins of angles
623 probRef = (1. - wAngle) * m_vecAngle_Time[ithAngle][ithJitter] +
624 wAngle * m_vecAngle_Time[jthAngle][ithJitter];
625 }
626 return jitter;
627}
628//+++++++++++++++++++++++++++++++++++++++++++++++
629bool TgcDigitMaker::efficiencyCheck(const TgcSensor sensor,
630 CLHEP::HepRandomEngine* rndmEngine) const {
631 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) < m_efficiency[sensor])
632 return true;
633 ATH_MSG_DEBUG("efficiencyCheck(): Hit removed for sensor= "
634 << sensor << "(0=WIRE,1=STRIP)");
635 return false;
636}
637//+++++++++++++++++++++++++++++++++++++++++++++++
638bool TgcDigitMaker::efficiencyCheck(const std::string& stationName,
639 const int stationEta, const int stationPhi,
640 const int gasGap, const TgcSensor sensor,
641 const double energyDeposit) const {
642 // If the energy deposit is equal to or greater than the threshold value of
643 // the chamber, return true.
644 return (energyDeposit >= getEnergyThreshold(stationName, stationEta,
645 stationPhi, gasGap, sensor));
646}
647//+++++++++++++++++++++++++++++++++++++++++++++++
648uint16_t TgcDigitMaker::bcTagging(const double digitTime, const double window,
649 const double offset) const {
650 const double calc_coll_time =
651 digitTime - offset; // calculated collision time
652 uint16_t bctag = 0;
653 if (-m_bunchCrossingTime < calc_coll_time &&
654 calc_coll_time < window - m_bunchCrossingTime) {
655 bctag |= 0x1;
656 }
657 if (0. < calc_coll_time && calc_coll_time < window) {
658 bctag |= 0x2;
659 }
660 if (m_bunchCrossingTime < calc_coll_time &&
661 calc_coll_time < window + m_bunchCrossingTime) {
662 bctag |= 0x4;
663 }
664 if (2. * m_bunchCrossingTime < calc_coll_time &&
665 calc_coll_time < window + 2. * m_bunchCrossingTime &&
666 m_doFourBunchDigitization) {
667 bctag |= 0x8;
668 }
669 return bctag;
670}
671//+++++++++++++++++++++++++++++++++++++++++++++++
672void TgcDigitMaker::addDigit(const Identifier id, const uint16_t bctag,
673 TgcDigitCollection* digits) {
674 for (int bc = TgcDigit::BC_PREVIOUS; bc <= TgcDigit::BC_NEXTNEXT; bc++) {
675 if ((bctag >> (bc - 1)) & 0x1) {
676 bool duplicate = false;
677 for (const auto digit : *digits) {
678 if (id == digit->identify() && digit->bcTag() == bc) {
679 duplicate = true;
680 break;
681 }
682 }
683 if (!duplicate) {
684 std::unique_ptr<TgcDigit> multihitDigit =
685 std::make_unique<TgcDigit>(id, bc);
686 digits->push_back(multihitDigit.release());
687 }
688 }
689 }
690 }
691
692StatusCode TgcDigitMaker::readFileOfEnergyThreshold() {
693 // Indices to be used
694 int iStationName, stationEta, stationPhi, gasGap, isStrip;
695
696 for (iStationName = 0; iStationName < N_STATIONNAME; iStationName++) {
697 for (stationEta = 0; stationEta < N_STATIONETA; stationEta++) {
698 for (stationPhi = 0; stationPhi < N_STATIONPHI; stationPhi++) {
699 for (gasGap = 0; gasGap < N_GASGAP; gasGap++) {
700 for (isStrip = 0; isStrip < N_ISSTRIP; isStrip++) {
701 m_energyThreshold[iStationName][stationEta][stationPhi]
702 [gasGap][isStrip] = -999999.;
703 }
704 }
705 }
706 }
707 }
708
709 // Find path to the TGC_Digitization_energyThreshold.dat file
710 const std::string fileName = "TGC_Digitization_energyThreshold.dat";
711 std::string fileWithPath = PathResolver::find_file(fileName, "DATAPATH");
712 if (fileWithPath.empty()) {
713 ATH_MSG_FATAL("readFileOfEnergyThreshold(): Could not find file "
714 << fileName);
715 return StatusCode::FAILURE;
716 }
717
718 // Open the TGC_Digitization_energyThreshold.dat file
719 std::ifstream ifs;
720 ifs.open(fileWithPath.c_str(), std::ios::in);
721 if (ifs.bad()) {
722 ATH_MSG_FATAL("readFileOfEnergyThreshold(): Could not open file "
723 << fileName);
724 return StatusCode::FAILURE;
725 }
726
727 double energyThreshold;
728 // Read the TGC_Digitization_energyThreshold.dat file
729 while (ifs.good()) {
730 ifs >> iStationName >> stationEta >> stationPhi >> gasGap >> isStrip >>
731 energyThreshold;
732 ATH_MSG_DEBUG("readFileOfEnergyThreshold"
733 << " stationName= " << iStationName << " stationEta= "
734 << stationEta << " stationPhi= " << stationPhi
735 << " gasGap= " << gasGap << " isStrip= " << isStrip
736 << " energyThreshold(MeV)= " << energyThreshold);
737
738 // Subtract offsets to use indices of energyThreshold array
739 iStationName -= OFFSET_STATIONNAME;
740 stationEta -= OFFSET_STATIONETA;
741 stationPhi -= OFFSET_STATIONPHI;
742 gasGap -= OFFSET_GASGAP;
743 isStrip -= OFFSET_ISSTRIP;
744
745 // Check the indices are valid
746 if (iStationName < 0 || iStationName >= N_STATIONNAME)
747 continue;
748 if (stationEta < 0 || stationEta >= N_STATIONETA)
749 continue;
750 if (stationPhi < 0 || stationPhi >= N_STATIONPHI)
751 continue;
752 if (gasGap < 0 || gasGap >= N_GASGAP)
753 continue;
754 if (isStrip < 0 || isStrip >= N_ISSTRIP)
755 continue;
756
757 m_energyThreshold[iStationName][stationEta][stationPhi][gasGap]
758 [isStrip] = energyThreshold;
759
760 // If it is the end of the file, get out from while loop.
761 if (ifs.eof())
762 break;
763 }
764
765 // Close the TGC_Digitization_energyThreshold.dat file
766 ifs.close();
767
768 return StatusCode::SUCCESS;
769}
770
771//--------------------------------------------
772unsigned int TgcDigitMaker::getRunPeriod() const
773{
774 // Used to determine the version of the TGC Dead Chambers text file to read.
775 // TODO There must be a better way of doing this -> ConditionsAlg?
776 SmartIF<IGeoModelSvc> geoModel{Gaudi::svcLocator()->service("GeoModelSvc")};
777 if (!geoModel) {
778 ATH_MSG_ERROR("getRunPeriod() Failed to find GeoModelSvc");
779 return 0;
780 }
781 std::string atlasVersion = geoModel->atlasVersion();
782
783 SmartIF<IRDBAccessSvc> rdbAccess{Gaudi::svcLocator()->service("RDBAccessSvc")};
784 if (!rdbAccess) {
785 ATH_MSG_ERROR("getRunPeriod() Failed to find RDBAccessSvc");
786 return 0;
787 }
788
789 IRDBRecordset_ptr atlasCommonRec =
790 rdbAccess->getRecordsetPtr("AtlasCommon", atlasVersion, "ATLAS");
791 unsigned int runperiod = 1;
792 if (atlasCommonRec->size() == 0)
793 runperiod = 1;
794 else {
795 std::string configVal = (*atlasCommonRec)[0]->getString("CONFIG");
796 if (configVal == "RUN1")
797 runperiod = 1;
798 else if (configVal == "RUN2")
799 runperiod = 2;
800 else if (configVal == "RUN3")
801 runperiod =
802 3; // currently runperiod 3 means no masking => ok for upgrade
803 else if (configVal == "RUN4")
804 runperiod =
805 3; // currently runperiod 3 means no masking => ok for upgrade
806 else {
807 runperiod = 0;
808 ATH_MSG_ERROR(
809 "Unexpected value for geometry config read from the database: "
810 << configVal);
811 }
812 }
813 return runperiod;
814}
815
816StatusCode TgcDigitMaker::readFileOfDeadChamber() {
817 // TODO There must be a better way of doing this -> ConditionsAlg?
818 // Indices to be used
819 int iStationName, stationEta, stationPhi, gasGap;
820
821 for (iStationName = 0; iStationName < N_STATIONNAME; iStationName++) {
822 for (stationEta = 0; stationEta < N_STATIONETA; stationEta++) {
823 for (stationPhi = 0; stationPhi < N_STATIONPHI; stationPhi++) {
824 for (gasGap = 0; gasGap < N_GASGAP; gasGap++) {
825 m_isDeadChamber[iStationName][stationEta][stationPhi]
826 [gasGap] = false;
827 }
828 }
829 }
830 }
831
832 unsigned int runperiod = getRunPeriod();
833 if (runperiod==0) {
834 ATH_MSG_FATAL("Could not determine run period.");
835 return StatusCode::FAILURE;
836 }
837
838 // Find path to the TGC_Digitization_deadChamber.dat file
839 std::string fileName;
840 if (runperiod == 1)
841 fileName = "TGC_Digitization_deadChamber.dat";
842 else if (runperiod == 2)
843 fileName = "TGC_Digitization_2016deadChamber.dat";
844 else if (runperiod == 3)
845 fileName = "TGC_Digitization_NOdeadChamber.dat";
846 else {
847 ATH_MSG_ERROR("Run Period " << runperiod
848 << " is unexpected in TgcDigitMaker - "
849 "using NOdeadChamber configuration.");
850 return StatusCode::FAILURE;
851 }
852 std::string fileWithPath = PathResolver::find_file(fileName, "DATAPATH");
853 if (fileWithPath.empty()) {
854 ATH_MSG_FATAL("readFileOfDeadChamber(): Could not find file "
855 << fileName);
856 return StatusCode::FAILURE;
857 }
858
859 // Open the TGC_Digitization_deadChamber.dat file
860 std::ifstream ifs;
861 ifs.open(fileWithPath.c_str(), std::ios::in);
862 if (ifs.bad()) {
863 ATH_MSG_FATAL("readFileOfDeadChamber(): Could not open file "
864 << fileName);
865 return StatusCode::FAILURE;
866 }
867
868 // Read the TGC_Digitization_deadChamber.dat file
869 unsigned int nDeadChambers = 0;
870 std::string comment;
871 while (ifs.good()) {
872 ifs >> iStationName >> stationEta >> stationPhi >> gasGap;
873 bool valid = getline(ifs, comment).good();
874 if (!valid)
875 break;
876
877 ATH_MSG_DEBUG("TgcDigitMaker::readFileOfDeadChamber"
878 << " stationName= " << iStationName << " stationEta= "
879 << stationEta << " stationPhi= " << stationPhi
880 << " gasGap= " << gasGap << " comment= " << comment);
881
882 // Subtract offsets to use indices of isDeadChamber array
883 iStationName -= OFFSET_STATIONNAME;
884 stationEta -= OFFSET_STATIONETA;
885 stationPhi -= OFFSET_STATIONPHI;
886 gasGap -= OFFSET_GASGAP;
887
888 // Check the indices are valid
889 if (iStationName < 0 || iStationName >= N_STATIONNAME)
890 continue;
891 if (stationEta < 0 || stationEta >= N_STATIONETA)
892 continue;
893 if (stationPhi < 0 || stationPhi >= N_STATIONPHI)
894 continue;
895 if (gasGap < 0 || gasGap >= N_GASGAP)
896 continue;
897
898 m_isDeadChamber[iStationName][stationEta][stationPhi][gasGap] = true;
899 nDeadChambers++;
900
901 // If it is the end of the file, get out from while loop.
902 if (ifs.eof())
903 break;
904 }
905
906 // Close the TGC_Digitization_deadChamber.dat file
907 ifs.close();
908
909 ATH_MSG_INFO("readFileOfDeadChamber: the number of dead chambers = "
910 << nDeadChambers);
911
912 return StatusCode::SUCCESS;
913}
914
915StatusCode TgcDigitMaker::readFileOfStripPosition() {
916 // Indices to be used
917 int iStationName, stationEta, channel;
918
919 for (iStationName = 0; iStationName < N_STATIONNAME; iStationName++) {
920 for (stationEta = 0; stationEta < N_ABSSTATIONETA; stationEta++) {
921 for (channel = 0; channel < N_STRIPCHANNEL; channel++) {
922 m_StripPos[iStationName][stationEta][channel] = 0.;
923 }
924 }
925 }
926
927 // Find path to the TGC_Digitization_StripPosition.dat file
928 const std::string fileName = "TGC_Digitization_StripPosition.dat";
929 std::string fileWithPath = PathResolver::find_file(fileName, "DATAPATH");
930 if (fileWithPath.empty()) {
931 ATH_MSG_FATAL("readFileOfStripPosition(): Could not find file "
932 << fileName);
933 return StatusCode::FAILURE;
934 }
935
936 // Open the TGC_Digitization_StripPosition.dat file
937 std::ifstream ifs;
938 ifs.open(fileWithPath.c_str(), std::ios::in);
939 if (ifs.bad()) {
940 ATH_MSG_FATAL("readFileOfStripPosition(): Could not open file "
941 << fileName);
942 return StatusCode::FAILURE;
943 }
944
945 // Read the TGC_Digitization_StripPosition.dat file
946 double strippos;
947 while (ifs.good()) {
948 ifs >> iStationName >> stationEta >> channel >> strippos;
949 ATH_MSG_DEBUG("readFileOfStripPosition"
950 << " stationName= " << iStationName << " stationEta= "
951 << stationEta << " channel= " << channel
952 << " StripPosition= " << strippos);
953
954 iStationName -= OFFSET_STATIONNAME;
955 stationEta -= OFFSET_ABSSTATIONETA;
956 channel -= OFFSET_STRIPCHANNEL;
957
958 // Check the indices are valid
959 if (iStationName < 0 || iStationName >= N_STATIONNAME)
960 continue;
961 if (stationEta < 0 || stationEta >= N_ABSSTATIONETA)
962 continue;
963 if (channel < 0 || channel >= N_STRIPCHANNEL)
964 continue;
965
966 m_StripPos[iStationName][stationEta][channel] = strippos;
967 // If it is the end of the file, get out from while loop.
968 if (ifs.eof())
969 break;
970 }
971 // Close the TGC_Digitization_StripPosition.dat file
972 ifs.close();
973
974 return StatusCode::SUCCESS;
975}
976
977double TgcDigitMaker::getEnergyThreshold(const std::string& stationName,
978 int stationEta, int stationPhi,
979 int gasGap,
980 const TgcSensor sensor) const {
981 // Convert std::string stationName to int iStationName from 41 to 48
982 int iStationName = getIStationName(stationName);
983
984 // Subtract offsets to use these as the indices of the energyThreshold array
985 iStationName -= OFFSET_STATIONNAME;
986 stationEta -= OFFSET_STATIONETA;
987 stationPhi -= OFFSET_STATIONPHI;
988 gasGap -= OFFSET_GASGAP;
989
990 double energyThreshold = +999999.;
991
992 // If the indices are valid, the energyThreshold array is fetched.
993 if ((iStationName >= 0 && iStationName < N_STATIONNAME) &&
994 (stationEta >= 0 && stationEta < N_STATIONETA) &&
995 (stationPhi >= 0 && stationPhi < N_STATIONPHI) &&
996 (gasGap >= 0 && gasGap < N_GASGAP)) {
997 energyThreshold = m_energyThreshold[iStationName][stationEta]
998 [stationPhi][gasGap][sensor];
999 }
1000
1001 // Show the energy threshold value
1002 ATH_MSG_VERBOSE("getEnergyThreshold"
1003 << " stationName= " << iStationName + OFFSET_STATIONNAME
1004 << " stationEta= " << stationEta + OFFSET_STATIONETA
1005 << " stationPhi= " << stationPhi + OFFSET_STATIONPHI
1006 << " gasGap= " << gasGap + OFFSET_GASGAP << " sensor= "
1007 << sensor << " energyThreshold(MeV)= " << energyThreshold);
1008
1009 return energyThreshold;
1010}
1011
1012//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1013void TgcDigitMaker::randomCrossTalk(
1014 const TgcDigitCrosstalkData* crosstalk, const Identifier elemId,
1015 const int gasGap, const TgcSensor sensor, const int channel,
1016 const float posInChan, const float digitTime, const float time_offset,
1017 CLHEP::HepRandomEngine* rndmEngine, TgcDigitCollection* digits) const {
1018 uint16_t station_number = m_idHelper->stationName(elemId);
1019 uint16_t station_eta = std::abs(m_idHelper->stationEta(elemId));
1020 uint16_t layer = gasGap;
1021 uint16_t layer_id = (station_number << 5) + (station_eta << 2) + layer;
1022
1023 if (station_number < OFFSET_STATIONNAME ||
1024 station_number >= OFFSET_STATIONNAME + N_STATIONNAME ||
1025 station_eta <= 0 || station_eta > N_ABSSTATIONETA || layer <= 0 ||
1026 layer > N_GASGAP) {
1027 ATH_MSG_ERROR("Unexpected indices are provided!");
1028 return;
1029 }
1030
1031 float prob1CrossTalk =
1032 this->getCrosstalkProbability(crosstalk, layer_id, sensor, 0);
1033 float prob11CrossTalk =
1034 this->getCrosstalkProbability(crosstalk, layer_id, sensor, 1);
1035 float prob20CrossTalk =
1036 this->getCrosstalkProbability(crosstalk, layer_id, sensor, 2);
1037 float prob21CrossTalk =
1038 this->getCrosstalkProbability(crosstalk, layer_id, sensor, 3);
1039
1040 int nCrossTalks_neg = 0;
1041 int nCrossTalks_pos = 0;
1042
1043 if (posInChan < prob1CrossTalk) {
1044 nCrossTalks_neg = 1; // 1-0
1045 } else if (posInChan > 1. - prob1CrossTalk) {
1046 nCrossTalks_pos = 1; // 0-1
1047 } else {
1048 double prob = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0);
1049 if (prob < prob11CrossTalk / (1. - 2. * prob1CrossTalk)) {
1050 nCrossTalks_neg = 1;
1051 nCrossTalks_pos = 1; // 1-1
1052 } else if (prob < (prob20CrossTalk + prob11CrossTalk) /
1053 (1. - 2. * prob1CrossTalk)) {
1054 if (posInChan < 0.5) {
1055 nCrossTalks_neg = 2;
1056 } // 2-0
1057 else {
1058 nCrossTalks_pos = 2;
1059 } // 0-2
1060 } else {
1061 if (prob <
1062 (prob20CrossTalk + prob11CrossTalk + 2. * prob21CrossTalk) /
1063 (1. - 2. * prob1CrossTalk)) {
1064 if (posInChan < 0.5) {
1065 nCrossTalks_neg = 2;
1066 nCrossTalks_pos = 1;
1067 } // 2-1
1068 else {
1069 nCrossTalks_neg = 1;
1070 nCrossTalks_pos = 2;
1071 } // 1-2
1072 }
1073 }
1074 }
1075
1076 if (nCrossTalks_neg == 0 && nCrossTalks_pos == 0)
1077 return; // No cross-talk case
1078
1079 // No time structure is implemented yet.
1080 float dt = digitTime;
1081 TgcStation station = (station_number > 46) ? kINNER : kOUTER;
1082 uint16_t bctag =
1083 bcTagging(dt, m_gateTimeWindow[station][sensor], time_offset);
1084 // obtain max channel number
1085 Identifier thisId =
1086 m_idHelper->channelID(elemId, gasGap, (int)sensor, channel);
1087 int maxChannelNumber = m_idHelper->channelMax(thisId);
1088
1089 for (int jChan = channel - nCrossTalks_neg;
1090 jChan <= channel + nCrossTalks_pos; jChan++) {
1091 if (jChan == channel || jChan < 1 || jChan > maxChannelNumber)
1092 continue;
1093
1094 Identifier newId =
1095 m_idHelper->channelID(elemId, gasGap, (int)sensor, jChan);
1096 addDigit(newId, bctag, digits); // TgcDigit can be duplicated.
1097 }
1098}
1099
1100bool TgcDigitMaker::isDeadChamber(const std::string& stationName,
1101 int stationEta, int stationPhi, int gasGap) {
1102 bool v_isDeadChamber = true;
1103
1104 // Convert std::string stationName to int iStationName from 41 to 48
1105 int iStationName = getIStationName(stationName);
1106
1107 // Subtract offsets to use these as the indices of the energyThreshold array
1108 iStationName -= OFFSET_STATIONNAME;
1109 stationEta -= OFFSET_STATIONETA;
1110 stationPhi -= OFFSET_STATIONPHI;
1111 gasGap -= OFFSET_GASGAP;
1112
1113 // If the indices are valid, the energyThreshold array is fetched.
1114 if ((iStationName >= 0 && iStationName < N_STATIONNAME) &&
1115 (stationEta >= 0 && stationEta < N_STATIONETA) &&
1116 (stationPhi >= 0 && stationPhi < N_STATIONPHI) &&
1117 (gasGap >= 0 && gasGap < N_GASGAP)) {
1118 v_isDeadChamber =
1119 m_isDeadChamber[iStationName][stationEta][stationPhi][gasGap];
1120 }
1121
1122 // Show the energy threshold value
1123 ATH_MSG_VERBOSE("TgcDigitMaker::getEnergyThreshold"
1124 << " stationName= " << iStationName + OFFSET_STATIONNAME
1125 << " stationEta= " << stationEta + OFFSET_STATIONETA
1126 << " stationPhi= " << stationPhi + OFFSET_STATIONPHI
1127 << " gasGap= " << gasGap + OFFSET_GASGAP
1128 << " isDeadChamber= " << v_isDeadChamber);
1129
1130 return v_isDeadChamber;
1131}
1132
1133int TgcDigitMaker::getIStationName(const std::string& stationName) {
1134 int iStationName = 0;
1135 if (stationName == "T1F")
1136 iStationName = 41;
1137 else if (stationName == "T1E")
1138 iStationName = 42;
1139 else if (stationName == "T2F")
1140 iStationName = 43;
1141 else if (stationName == "T2E")
1142 iStationName = 44;
1143 else if (stationName == "T3F")
1144 iStationName = 45;
1145 else if (stationName == "T3E")
1146 iStationName = 46;
1147 else if (stationName == "T4F")
1148 iStationName = 47;
1149 else if (stationName == "T4E")
1150 iStationName = 48;
1151
1152 return iStationName;
1153}
1154
1155float TgcDigitMaker::getStripPosition(const std::string& stationName,
1156 int stationEta, int channel) const {
1157 // Convert std::string stationName to int iStationName from 41 to 48
1158 int iStationName = getIStationName(stationName);
1159
1160 // Subtract offsets to use these as the indices of the energyThreshold array
1161 iStationName -= OFFSET_STATIONNAME;
1162 stationEta -= OFFSET_ABSSTATIONETA;
1163 channel -= OFFSET_STRIPCHANNEL;
1164
1165 // Check the indices are valid
1166 if (iStationName < 0 || iStationName >= N_STATIONNAME)
1167 return 0.;
1168 if (stationEta < 0 || stationEta >= N_ABSSTATIONETA)
1169 return 0.;
1170 if (channel < 0 || channel >= N_STRIPCHANNEL)
1171 return 0.;
1172
1173 return m_StripPos[iStationName][stationEta][channel];
1174}
1175
1176float TgcDigitMaker::timeDiffByCableRadiusOfInner(const int iStationName,
1177 const int stationPhi,
1178 const int channel) {
1179 float delay{0.};
1180 if (iStationName != 47 && iStationName != 48)
1181 return delay; // Big Wheel has no delay for this.
1182
1183 if (channel < 12 || (channel > 16 && channel < 27)) {
1184 int cablenum = (stationPhi >= 13) ? 25 - stationPhi : stationPhi;
1185 delay += 2.3 * CLHEP::ns - 0.06 * CLHEP::ns * float(cablenum);
1186 }
1187 return delay;
1188}
1189
1190double TgcDigitMaker::getSigPropTimeDelay(const float cableDistance) {
1191 constexpr std::array<double, 2> par{0.0049 * CLHEP::ns, 0.0002 * CLHEP::ns};
1192 return par[0] * std::pow(cableDistance / CLHEP::m, 2) +
1193 par[1] * cableDistance / CLHEP::m;
1194}
1195
1196float TgcDigitMaker::getDistanceToAsdFromSensor(
1197 const TgcDigitASDposData* readCdo, const int iStationName,
1198 const int stationEta, const int stationPhi, const TgcSensor sensor,
1199 const int channel, const float position) const {
1200 // EIFI has different parameter for phi, BW is same parameter for phi (i.e.
1201 // -99 in DB).
1202 int phiId = (iStationName >= 47) ? stationPhi : 0x1f;
1203 uint16_t chamberId = (iStationName << 8) + (stationEta << 5) + phiId;
1204
1205 unsigned int asdNo = static_cast<unsigned int>(channel - 1) /
1206 TgcDigitASDposData::N_CHANNELINPUT_TOASD;
1207
1208 float asd_position = 0.;
1209 const auto& map = (sensor == kSTRIP) ? readCdo->stripAsdPos : readCdo->wireAsdPos;
1210 auto it = map.find(chamberId);
1211
1212 if (it != map.end()) {
1213 asd_position = it->second[asdNo] * CLHEP::m;
1214 } else {
1215 ATH_MSG_WARNING("Undefined chamberID is provided! station="
1216 << iStationName << ", eta=" << stationEta
1217 << ", phi=" << phiId);
1218 }
1219
1220 float distance = fabs(position - asd_position);
1221 return distance;
1222}
1223
1224float TgcDigitMaker::getTimeOffset(const TgcDigitTimeOffsetData* readCdo,
1225 const uint16_t station_num,
1226 const int station_eta,
1227 const TgcSensor sensor) {
1228 uint16_t chamberId =
1229 (station_num << 3) + static_cast<uint16_t>(std::abs(station_eta));
1230 return ((sensor == TgcSensor::kSTRIP)
1231 ? readCdo->stripOffset.find(chamberId)->second
1232 : readCdo->wireOffset.find(chamberId)->second);
1233}
1234
1235float TgcDigitMaker::getCrosstalkProbability(
1236 const TgcDigitCrosstalkData* readCdo, const uint16_t layer_id,
1237 const TgcSensor sensor, const unsigned int index_prob) {
1238 if (readCdo == nullptr)
1239 return 0.; // no crosstalk
1240 return ((sensor == TgcSensor::kSTRIP)
1241 ? readCdo->getStripProbability(layer_id, index_prob)
1242 : readCdo->getWireProbability(layer_id, index_prob));
1243}
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
IRDBAccessSvc.h
Definition of the abstract IRDBAccessSvc interface.
IRDBRecordset_ptr
std::shared_ptr< IRDBRecordset > IRDBRecordset_ptr
Definition IRDBAccessSvc.h:27
IRDBRecord.h
Definition of the abstract IRDBRecord interface.
IRDBRecordset.h
Definition of the abstract IRDBRecordset interface.
delay
double delay(std::size_t d)
Definition JetTrigTimerTest.cxx:11
bins
static const std::vector< std::string > bins
Definition MuonTriggerSFFilesTest.cxx:22
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition TRTDetectorFactory_Full.cxx:71
AthMessaging::msg
MsgStream & msg() const
The standard message stream.
Definition AthMessaging.h:167
AthMessaging::msgLvl
bool msgLvl(const MSG::Level lvl) const
Test the output level.
Definition AthMessaging.h:151
AthMessaging::AthMessaging
AthMessaging(IMessageSvc *msgSvc, const std::string &name)
Constructor.
Definition AthMessaging.cxx:13
IRDBRecordset::size
virtual unsigned int size() const =0
IdContext
This class saves the "context" of an expanded identifier (ExpandedIdentifier) for compact or hash ver...
Definition IdContext.h:26
IdContext::begin_index
size_type begin_index() const
Definition IdContext.h:45
IdContext::end_index
size_type end_index() const
Definition IdContext.h:46
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
MuonGM::MuonDetectorManager
The MuonDetectorManager stores the transient representation of the Muon Spectrometer geometry and pro...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MuonDetectorManager.h:51
MuonGM::TgcReadoutElement
A TgcReadoutElement corresponds to a single TGC chamber; therefore typically a TGC station contains s...
Definition MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/TgcReadoutElement.h:42
MuonGM::TgcReadoutElement::chamberWidth
double chamberWidth(double z) const
MuonGM::TgcReadoutElement::nWires
int nWires(int gasGap) const
Returns the total number of wires in a given gas gap.
MuonGM::TgcReadoutElement::nStrips
int nStrips(int gasGap) const
Returns the number of strips in a given gas gap.
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path)
Definition PathResolver.cxx:220
TGCSimHit::TGCid
int TGCid() const
Definition TGCSimHit.h:45
TGCSimHit::energyDeposit
double energyDeposit() const
Definition TGCSimHit.h:46
TGCSimHit::localDireCos
const Amg::Vector3D & localDireCos() const
Definition TGCSimHit.h:44
TGCSimHit::localPosition
const Amg::Vector3D & localPosition() const
Definition TGCSimHit.h:43
TgcDigitASDposData::stripAsdPos
std::map< uint16_t, std::vector< float > > stripAsdPos
Definition TgcDigitASDposData.h:35
TgcDigitASDposData::wireAsdPos
std::map< uint16_t, std::vector< float > > wireAsdPos
Definition TgcDigitASDposData.h:36
TgcDigitCrosstalkData::getStripProbability
float getStripProbability(const uint16_t layer_id, const unsigned int index_prob) const
Definition TgcDigitCrosstalkData.cxx:15
TgcDigitCrosstalkData::getWireProbability
float getWireProbability(const uint16_t layer_id, const unsigned int index_prob) const
Definition TgcDigitCrosstalkData.cxx:28
TgcDigitMaker::getTimeOffset
static float getTimeOffset(const TgcDigitTimeOffsetData *readCdo, const uint16_t station_num, const int station_eta, const TgcSensor sensor)
Method to get time offset to absorb signal delay.
Definition TgcDigitMaker.cxx:1224
TgcDigitMaker::getCrosstalkProbability
static float getCrosstalkProbability(const TgcDigitCrosstalkData *readCdo, const uint16_t layer_id, const TgcSensor sensor, const unsigned int index_prob)
Method to get the channel crosstalk probability.
Definition TgcDigitMaker.cxx:1235
TgcDigitMaker::readFileOfDeadChamber
StatusCode readFileOfDeadChamber()
Read share/TGC_Digitization_deadChamber.dat file.
Definition TgcDigitMaker.cxx:816
TgcDigitMaker::initialize
StatusCode initialize()
Initializes TgcHitIdHelper, TgcIdHelper and random number of a stream for the digitization.
Definition TgcDigitMaker.cxx:59
TgcDigitMaker::getSigPropTimeDelay
static double getSigPropTimeDelay(const float cableDistance)
Method to get signal propagation time delay.
Definition TgcDigitMaker.cxx:1190
TgcDigitMaker::addDigit
static void addDigit(const Identifier id, const uint16_t bctag, TgcDigitCollection *digits)
Definition TgcDigitMaker.cxx:672
TgcDigitMaker::efficiencyCheck
bool efficiencyCheck(const TgcSensor sensor, CLHEP::HepRandomEngine *rndmEngine) const
Determines whether a hit is detected or not.
Definition TgcDigitMaker.cxx:629
TgcDigitMaker::m_efficiency
float m_efficiency[N_SENSOR]
Definition TgcDigitMaker.h:196
TgcDigitMaker::readFileOfTimeJitter
StatusCode readFileOfTimeJitter()
Reads parameters for intrinsic time response from timejitter.dat.
Definition TgcDigitMaker.cxx:547
TgcDigitMaker::getDistanceToAsdFromSensor
float getDistanceToAsdFromSensor(const TgcDigitASDposData *readCdo, const int iStationName, const int stationEta, const int stationPhi, const TgcSensor sensor, const int channel, const float position) const
Method to get propagation time to the ASD from the sensor.
Definition TgcDigitMaker.cxx:1196
TgcDigitMaker::getIStationName
static int getIStationName(const std::string &staionName)
Get stationName integer from stationName string.
Definition TgcDigitMaker.cxx:1133
TgcDigitMaker::timeDiffByCableRadiusOfInner
static float timeDiffByCableRadiusOfInner(const int iStationName, const int stationPhi, const int channel)
Method to get time difference by cable radius of inner.
Definition TgcDigitMaker.cxx:1176
TgcDigitMaker::getStripPosition
float getStripPosition(const std::string &stationName, int stationEta, int channel) const
Method to get position of Strip channel.
Definition TgcDigitMaker.cxx:1155
TgcDigitMaker::m_bunchCrossingTime
double m_bunchCrossingTime
Definition TgcDigitMaker.h:205
TgcDigitMaker::m_vecAngle_Time
std::vector< std::vector< float > > m_vecAngle_Time
Definition TgcDigitMaker.h:190
TgcDigitMaker::getRunPeriod
unsigned int getRunPeriod() const
Determine the run period.
Definition TgcDigitMaker.cxx:772
TgcDigitMaker::isDeadChamber
bool isDeadChamber(const std::string &stationName, int stationEta, int stationPhi, int gasGap)
Method to check a chamber is dead or active.
Definition TgcDigitMaker.cxx:1100
TgcDigitMaker::timeJitter
float timeJitter(const Amg::Vector3D &, CLHEP::HepRandomEngine *rndmEngine) const
Calculates intrinsic time response according to incident angle of a track based on time response para...
Definition TgcDigitMaker.cxx:597
TgcDigitMaker::readFileOfEnergyThreshold
StatusCode readFileOfEnergyThreshold()
Read share/TGC_Digitization_energyThreshold.dat file.
Definition TgcDigitMaker.cxx:692
TgcDigitMaker::m_energyThreshold
double m_energyThreshold[N_STATIONNAME][N_STATIONETA][N_STATIONPHI][N_GASGAP][N_ISSTRIP]
Energy threshold value for each chamber.
Definition TgcDigitMaker.h:183
TgcDigitMaker::randomCrossTalk
void randomCrossTalk(const TgcDigitCrosstalkData *crosstalk, const Identifier elemId, const int gasGap, const TgcSensor sensor, const int channel, const float posInStrip, const float digitTime, const float time_offset, CLHEP::HepRandomEngine *rndmEngine, TgcDigitCollection *digits) const
Definition TgcDigitMaker.cxx:1013
TgcDigitMaker::readFileOfStripPosition
StatusCode readFileOfStripPosition()
Read share/TGC_Digitization_StripPosition.dat file.
Definition TgcDigitMaker.cxx:915
TgcDigitMaker::m_gateTimeWindow
double m_gateTimeWindow[N_STATION][N_SENSOR]
define the time windows for signals from wiregangs and strips.
Definition TgcDigitMaker.h:204
TgcDigitMaker::TgcDigitMaker
TgcDigitMaker(const TgcHitIdHelper *hitIdHelper, const MuonGM::MuonDetectorManager *mdManager, const bool doFourBunch)
Definition TgcDigitMaker.cxx:37
TgcDigitMaker::m_mdManager
const MuonGM::MuonDetectorManager * m_mdManager
Definition TgcDigitMaker.h:194
TgcDigitMaker::executeDigi
TgcDigitCollection * executeDigi(const TGCSimHit *hit, const double globalHitTime, const TgcDigitASDposData *ASDpos, const TgcDigitTimeOffsetData *TOffset, const TgcDigitCrosstalkData *Crosstalk, CLHEP::HepRandomEngine *rndmEngine)
A single hit can be digitized in the two directions independently: R and phi directions.
Definition TgcDigitMaker.cxx:88
TgcDigitMaker::m_doFourBunchDigitization
bool m_doFourBunchDigitization
Activate four bunch digitization.
Definition TgcDigitMaker.h:208
TgcDigitMaker::getEnergyThreshold
double getEnergyThreshold(const std::string &stationName, int stationEta, int stationPhi, int gasGap, const TgcSensor sensor) const
Get energy threshold value for each chamber.
Definition TgcDigitMaker.cxx:977
TgcDigitMaker::bcTagging
uint16_t bcTagging(const double digittime, const double window, const double offset) const
Definition TgcDigitMaker.cxx:648
TgcDigitMaker::m_isDeadChamber
bool m_isDeadChamber[N_STATIONNAME][N_STATIONETA][N_STATIONPHI][N_GASGAP]
Dead chamber flag for each chamber.
Definition TgcDigitMaker.h:185
TgcDigitMaker::m_hitIdHelper
const TgcHitIdHelper * m_hitIdHelper
Definition TgcDigitMaker.h:192
TgcDigitMaker::m_idHelper
const TgcIdHelper * m_idHelper
Definition TgcDigitMaker.h:195
TgcDigitMaker::m_StripPos
float m_StripPos[N_STATIONNAME][N_ABSSTATIONETA][N_STRIPCHANNEL]
Position of Strip Channel (Longer base or Shorter base)
Definition TgcDigitMaker.h:188
TgcDigitTimeOffsetData::stripOffset
std::map< uint16_t, float > stripOffset
Definition TgcDigitTimeOffsetData.h:25
TgcDigitTimeOffsetData::wireOffset
std::map< uint16_t, float > wireOffset
Definition TgcDigitTimeOffsetData.h:26
TgcDigit::BC_PREVIOUS
@ BC_PREVIOUS
Definition TgcDigit.h:37
TgcDigit::BC_NEXTNEXT
@ BC_NEXTNEXT
Definition TgcDigit.h:37
TgcHitIdHelper::GetHelper
static const TgcHitIdHelper * GetHelper()
Definition TgcHitIdHelper.cxx:23
map
STL class.
verbose
bool verbose
Definition hcg.cxx:73
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Generated on for ATLAS Offline Software by doxygen 1.14.0