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