13 #include "CLHEP/Random/RandFlat.h"
14 #include "CLHEP/Random/RandomEngine.h"
15 #include "CLHEP/Units/SystemOfUnits.h"
16 #include "GaudiKernel/MsgStream.h"
33 unsigned int runperiod,
const bool doFourBunch)
34 :
AthMessaging(
"TgcDigitMaker"), m_doFourBunchDigitization(doFourBunch) {
81 return StatusCode::SUCCESS;
88 const TGCSimHit* hit,
const double globalHitTime,
91 CLHEP::HepRandomEngine* rndmEngine) {
93 constexpr
float sensor_propagation_time =
101 constexpr
float wire_pitch = 1.8 *
CLHEP::mm;
102 constexpr
float zwidth_frame = 17. *
CLHEP::mm;
107 int Id = hit->
TGCid();
134 <<
"context begin_index = " << tgcContext.
begin_index()
135 <<
" context end_index = " << tgcContext.
end_index()
136 <<
" the identifier is " << elemId);
140 std::unique_ptr<TgcDigitCollection> digits =
141 std::make_unique<TgcDigitCollection>(elemId, coll_hash);
144 float height = tgcChamber->
getRsize();
145 float hmin = sqrt(
pow(centreChamber.x(), 2) +
pow(centreChamber.y(), 2)) -
155 float distanceZ = 1.4 *
CLHEP::mm / direCos[0] * direCos[2];
156 float zLocal = localPos.z() + distanceZ;
159 float distanceY = 1.4 *
CLHEP::mm / direCos[0] * direCos[1];
163 float yLocal = ySign * (localPos.y() + distanceY);
168 double tofCorrection =
169 (sqrt(
pow(hmin + zwidth_frame, 2) +
pow(centreChamber.z(), 2)) /
173 float bunchTime = globalHitTime - tofCorrection;
175 static const float jitterInitial = 9999.;
176 float jitter = jitterInitial;
197 kWIRE, rndmEngine))) {
200 float posInWireGroup[2] = {0., 0.};
201 for (
int iPosition = 0; iPosition < 2; iPosition++) {
202 int nWireOffset = (std::abs(
stationEta) == 5 ||
210 double zPosInSensArea =
211 zLocal +
static_cast<double>(tgcChamber->
nWires(ilyr) -
216 if (zPosInSensArea < 0. ||
217 zPosInSensArea > tgcChamber->
nWires(ilyr) * wire_pitch) {
218 iWireGroup[iPosition] = 0;
219 posInWireGroup[iPosition] = 0.;
221 "executeDigi(): Wire: Hit position located at outside of a "
225 <<
" position: " << zPosInSensArea <<
" xlocal: " << zLocal
226 <<
" dir cosine: " << direCos[0] <<
"/" << direCos[1] <<
"/"
228 <<
" active region: " << height - zwidth_frame * 2.);
232 while (wire_pitch * (
static_cast<float>(wire_index)) <
234 igang <= tgcChamber->nWireGangs(ilyr)) {
235 wire_index += tgcChamber->
nWires(ilyr, igang);
238 posInWireGroup[iPosition] =
239 (zPosInSensArea / wire_pitch -
240 (
static_cast<float>(wire_index))) /
242 tgcChamber->
nWires(ilyr, igang - 1))) +
245 iWireGroup[iPosition] = ((1 == igang) ? 1 : igang - 1);
249 unsigned int jWG[2] = {
250 (iWireGroup[0] <= iWireGroup[1]) ? (
unsigned int)(0)
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]);
262 for (
int iwg = iWG[0]; iwg <= iWG[1]; iwg++) {
263 if (1 <= iwg && iwg <= tgcChamber->nWireGangs(ilyr)) {
265 float wire_timeOffset =
271 if (iStationName > 46)
279 float ySignPhi = (
stationPhi % 2 == 1) ? -1. : +1.;
283 float wTimeDiffByRadiusOfInner =
287 bunchTime + jitter + wTimeDiffByRadiusOfInner;
288 float wASDDis{-1000.};
289 if (ASDpos !=
nullptr) {
294 wASDDis + (ySignPhi * yLocal +
297 sensor_propagation_time *
300 cable_propagation_time * wASDDis;
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 *
320 cable_propagation_time * wASDDis
321 <<
" time difference by cable radius of inner station: "
322 << wTimeDiffByRadiusOfInner);
327 addDigit(newId, bctag, digits.get());
331 posInWG[0], digit_time, wire_timeOffset,
332 rndmEngine, digits.get());
336 "WireGroup: newid breakdown digitTime x/y/z direcos "
337 "height_gang bctag: "
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) <<
" "
349 "Wire Gang id is out of range. id, x/y/z, height: "
350 << iwg <<
" " << localPos.x() <<
"/" << localPos.y() <<
"/"
351 << localPos.z() <<
" " << height);
359 if ((ilyr != 2 || (
stationName.compare(0, 2,
"T1") !=
363 sensor, rndmEngine)) ||
372 float posInStrip[2] = {0., 0.};
374 for (
int iPosition = 0; iPosition < 2; iPosition++) {
377 float zPos = zLocal + height / 2. - zwidth_frame;
379 iStrip[iPosition] = 0;
380 posInStrip[iPosition] = 0.;
382 "Strip: Hit position located at outside of a sensitive "
383 "volume, id, position, xlocal0/1, dir cosine: "
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.;
391 "Strip: Hit position located at outside of a sensitive "
392 "volume, id, position, active region: "
394 <<
"/" << ilyr << zPos <<
" "
395 << height - zwidth_frame * 2.);
404 float phiLocal = atan2(yLocal, zLocal + height / 2. + hmin);
407 "dphi, phiLocal, yLocal, zLocall+ height/2.+hmin: "
408 << dphi <<
" " << phiLocal <<
" " << yLocal <<
" "
409 << zLocal + height / 2. + hmin);
414 istr =
static_cast<int>(floor(phiLocal / dphi + 15.75)) + 1;
415 posInStrip[iPosition] =
416 phiLocal / dphi + 15.75 -
static_cast<float>(istr - 1);
418 istr =
static_cast<int>(floor(
419 (phiLocal - dphi * 14.25) / (dphi / 2.))) +
421 posInStrip[iPosition] =
422 (phiLocal - dphi * 14.25) / (dphi / 2.) -
423 static_cast<float>(istr - 31);
426 istr =
static_cast<int>(floor(phiLocal / dphi + 16.25)) + 1;
427 posInStrip[iPosition] =
428 phiLocal / dphi + 16.25 -
static_cast<float>(istr - 1);
430 istr =
static_cast<int>(floor(
431 (phiLocal + dphi * 14.25) / (dphi / 2.))) +
433 posInStrip[iPosition] =
434 (phiLocal + dphi * 14.25) / (dphi / 2.) -
435 static_cast<float>(istr - 3);
440 posInStrip[iPosition] = 0.;
441 }
else if (32 < istr) {
443 posInStrip[iPosition] = 0.;
445 iStrip[iPosition] = istr;
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]) {
460 float strip_timeOffset =
465 for (
int istr = iStr[0]; istr <= iStr[1]; istr++) {
466 if (1 <= istr && istr <= 32) {
468 if (jitter > jitterInitial - 0.1) {
477 float sTimeDiffByRadiusOfInner =
482 sensor_propagation_time *
483 (height / 2. + zwidth_frame + zSignEta * zLocal) +
484 sTimeDiffByRadiusOfInner;
485 float sASDDis{-1000};
486 if (ASDpos !=
nullptr) {
491 float sCableDis = sASDDis + (height / 2. + zwidth_frame +
494 sASDDis * cable_propagation_time;
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);
516 addDigit(newId, bctag, digits.get());
518 if (istr == iStr[0]) {
520 iStr[0], posInStr[0], sDigitTime,
521 strip_timeOffset, rndmEngine,
526 "Strip: newid breakdown digitTime x/y/z direcos "
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);
538 << ilyr <<
" " << istr);
543 return digits.release();
548 const char*
const fileName =
"TGC_Digitization_timejitter.dat";
552 if (!fileWithPath.empty()) {
553 ifs.open(fileWithPath.c_str(), std::ios::in);
557 return StatusCode::FAILURE;
563 return StatusCode::FAILURE;
577 "readFileOfTimeJitter(): Timejitter, angle, Number of bins, prob. "
581 for (
int j = 0; j < 41 ; j++) {
594 return StatusCode::SUCCESS;
598 CLHEP::HepRandomEngine* rndmEngine)
const {
599 float injectionAngle =
600 atan2(fabs(direCosLocal[2]), fabs(direCosLocal[0])) /
CLHEP::degree;
602 int ithAngle =
static_cast<int>(injectionAngle / 5.);
603 float wAngle = injectionAngle / 5. -
static_cast<float>(ithAngle);
609 jthAngle = ithAngle + 1;
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. *
620 int ithJitter =
static_cast<int>(jitter);
630 CLHEP::HepRandomEngine* rndmEngine)
const {
631 if (CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0) <
m_efficiency[sensor])
634 << sensor <<
"(0=WIRE,1=STRIP)");
649 const double offset)
const {
650 const double calc_coll_time =
657 if (0. < calc_coll_time && calc_coll_time < window) {
675 if ((bctag >> (bc - 1)) & 0
x1) {
676 bool duplicate =
false;
677 for (
const auto digit : *digits) {
678 if (
id ==
digit->identify() &&
digit->bcTag() == bc) {
684 std::unique_ptr<TgcDigit> multihitDigit =
685 std::make_unique<TgcDigit>(
id, bc);
686 digits->push_back(multihitDigit.release());
696 for (iStationName = 0; iStationName <
N_STATIONNAME; iStationName++) {
710 const std::string
fileName =
"TGC_Digitization_energyThreshold.dat";
712 if (fileWithPath.empty()) {
713 ATH_MSG_FATAL(
"readFileOfEnergyThreshold(): Could not find file "
715 return StatusCode::FAILURE;
720 ifs.open(fileWithPath.c_str(), std::ios::in);
722 ATH_MSG_FATAL(
"readFileOfEnergyThreshold(): Could not open file "
724 return StatusCode::FAILURE;
727 double energyThreshold;
733 <<
" stationName= " << iStationName <<
" stationEta= "
736 <<
" energyThreshold(MeV)= " << energyThreshold);
752 if (gasGap < 0 || gasGap >=
N_GASGAP)
768 return StatusCode::SUCCESS;
775 for (iStationName = 0; iStationName <
N_STATIONNAME; iStationName++) {
789 fileName =
"TGC_Digitization_deadChamber.dat";
791 fileName =
"TGC_Digitization_2016deadChamber.dat";
793 fileName =
"TGC_Digitization_NOdeadChamber.dat";
796 <<
" is unexpected in TgcDigitMaker - "
797 "using NOdeadChamber configuration.");
798 return StatusCode::FAILURE;
801 if (fileWithPath.empty()) {
802 ATH_MSG_FATAL(
"readFileOfDeadChamber(): Could not find file "
804 return StatusCode::FAILURE;
809 ifs.open(fileWithPath.c_str(), std::ios::in);
811 ATH_MSG_FATAL(
"readFileOfDeadChamber(): Could not open file "
813 return StatusCode::FAILURE;
817 unsigned int nDeadChambers = 0;
826 <<
" stationName= " << iStationName <<
" stationEta= "
843 if (gasGap < 0 || gasGap >=
N_GASGAP)
857 ATH_MSG_INFO(
"readFileOfDeadChamber: the number of dead chambers = "
860 return StatusCode::SUCCESS;
867 for (iStationName = 0; iStationName <
N_STATIONNAME; iStationName++) {
876 const std::string
fileName =
"TGC_Digitization_StripPosition.dat";
878 if (fileWithPath.empty()) {
879 ATH_MSG_FATAL(
"readFileOfStripPosition(): Could not find file "
881 return StatusCode::FAILURE;
886 ifs.open(fileWithPath.c_str(), std::ios::in);
888 ATH_MSG_FATAL(
"readFileOfStripPosition(): Could not open file "
890 return StatusCode::FAILURE;
898 <<
" stationName= " << iStationName <<
" stationEta= "
900 <<
" StripPosition= " << strippos);
922 return StatusCode::SUCCESS;
938 double energyThreshold = +999999.;
955 << sensor <<
" energyThreshold(MeV)= " << energyThreshold);
957 return energyThreshold;
964 const float posInChan,
const float digitTime,
const float time_offset,
979 float prob1CrossTalk =
981 float prob11CrossTalk =
983 float prob20CrossTalk =
985 float prob21CrossTalk =
988 int nCrossTalks_neg = 0;
989 int nCrossTalks_pos = 0;
991 if (posInChan < prob1CrossTalk) {
993 }
else if (posInChan > 1. - prob1CrossTalk) {
996 double prob = CLHEP::RandFlat::shoot(rndmEngine, 0.0, 1.0);
997 if (
prob < prob11CrossTalk / (1. - 2. * prob1CrossTalk)) {
1000 }
else if (
prob < (prob20CrossTalk + prob11CrossTalk) /
1001 (1. - 2. * prob1CrossTalk)) {
1002 if (posInChan < 0.5) {
1003 nCrossTalks_neg = 2;
1006 nCrossTalks_pos = 2;
1010 (prob20CrossTalk + prob11CrossTalk + 2. * prob21CrossTalk) /
1011 (1. - 2. * prob1CrossTalk)) {
1012 if (posInChan < 0.5) {
1013 nCrossTalks_neg = 2;
1014 nCrossTalks_pos = 1;
1017 nCrossTalks_neg = 1;
1018 nCrossTalks_pos = 2;
1024 if (nCrossTalks_neg == 0 && nCrossTalks_pos == 0)
1028 float dt = digitTime;
1037 for (
int jChan =
channel - nCrossTalks_neg;
1038 jChan <=
channel + nCrossTalks_pos; jChan++) {
1039 if (jChan ==
channel || jChan < 1 || jChan > maxChannelNumber)
1050 bool v_isDeadChamber =
true;
1076 <<
" isDeadChamber= " << v_isDeadChamber);
1078 return v_isDeadChamber;
1082 int iStationName = 0;
1100 return iStationName;
1128 if (iStationName != 47 && iStationName != 48)
1147 const int channel,
const float position)
const {
1150 int phiId = (iStationName >= 47) ?
stationPhi : 0x1f;
1153 unsigned int asdNo =
static_cast<unsigned int>(
channel - 1) /
1156 float asd_position = 0.;
1158 auto it = map.find(chamberId);
1160 if (
it != map.end()) {
1165 <<
", phi=" << phiId);
1168 float distance = fabs(position - asd_position);
1177 (station_num << 3) + static_cast<uint16_t>(std::abs(
station_eta));
1178 return ((sensor == TgcSensor::kSTRIP)
1180 : readCdo->
wireOffset.find(chamberId)->second);
1185 const TgcSensor sensor,
const unsigned int index_prob) {
1186 if (readCdo ==
nullptr)
1188 return ((sensor == TgcSensor::kSTRIP)