80 const Identifier layid = idHelper.channelID(offlineId,
81 idHelper.multilayer(offlineId),
82 idHelper.gasGap(offlineId),
83 sTgcIdHelper::sTgcChannelTypes::Wire, 1);
93 if( !detEl ){
return {}; }
99 constexpr
double length_tolerance = 0.01;
115 Amg::Vector2D posOnSurf_wire{0.5 * (hitOnSurface_wire + pre_pos_wire_surf).block<2,1>(0,0)};
132 if ((wire_number < 1) || (wire_number > number_wires)) {
136 if ((wire_number < 1) || (wire_number > number_wires)) {
142 Ionization ionization =
pointClosestApproach(detEl, layid, wire_number, pre_pos_wire_surf, hitOnSurface_wire);
143 double dist_wire = ionization.distance;
144 if (dist_wire > 0.) {
149 if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {
adjacent = -1;}
154 double dist_wire_adj = ion_adj.distance;
155 if ((dist_wire_adj > 0.) && (dist_wire_adj < dist_wire)) {
156 dist_wire = dist_wire_adj;
158 ionization = std::move(ion_adj);
161 ATH_MSG_DEBUG(
"Failed to get the distance between the wire number = " << wire_number
163 <<
". Number of wires = " << number_wires
169 posOnSurf_wire = ionization.posOnWire.block<2,1>(0,0);
174 <<
" posOnTrack: " <<
Amg::toString(ionization.posOnSegment, 3)
187 const double wire_pitch = detEl->
wirePitch();
188 if ((dist_wire > 0.) && (std::abs(hit.
particleEncoding()) == 13) && (dist_wire > (wire_pitch/2))) {
189 ATH_MSG_DEBUG(
"Distance to the nearest wire (" << dist_wire <<
") is greater than expected.");
199 if (dist_wire > wire_pitch) {
205 const double par_kappa = gamParam.kParameter;
206 const double par_theta = gamParam.thetaParameter;
209 double gamma_mpv = (par_kappa - 1) * par_theta;
211 if (gamma_mpv < 0.) {gamma_mpv = 0.;}
212 const double t0_par = most_prob_time - gamma_mpv;
219 double digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
225 constexpr
int shoot_limit = 4;
226 int shoot_counter = 0;
227 while (digit_time < 0.) {
228 if (shoot_counter > shoot_limit) {
232 digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
237 <<
", time = " << digit_time
238 <<
", k parameter = " << par_kappa
239 <<
", theta parameter = " << par_theta
240 <<
", most probable time = " << most_prob_time);
244 if (condContainers.efficiencies) {
245 const double efficiency = condContainers.efficiencies->getEfficiency(layid);
247 if (CLHEP::RandFlat::shoot(condContainers.rndmEngine,0.0,1.0) >
efficiency)
return {};
254 const double stripPropagationTime = 0.;
256 double sDigitTimeWire = digit_time;
257 double sDigitTimePad = sDigitTimeWire;
258 double sDigitTimeStrip = sDigitTimeWire + stripPropagationTime;
266 int channelType = sTgcIdHelper::sTgcChannelTypes::Strip;
268 Identifier newId = idHelper.channelID(layid,idHelper.multilayer(layid),
269 idHelper.gasGap(layid), channelType, 1,
isValid);
277 const Amg::Vector2D posOnSurf_strip(hitOnSurface_strip.x(),hitOnSurface_strip.y());
278 bool insideBounds = SURF_STRIP.
insideBounds(posOnSurf_strip);
286 int stripNumber = detEl->
stripNumber(posOnSurf_strip, newId);
287 if( stripNumber == -1 ){
289 const double new_posX = (posOnSurf_strip.x() > 0.0)? posOnSurf_strip.x() - length_tolerance
290 : posOnSurf_strip.x() + length_tolerance;
291 const Amg::Vector2D new_position(new_posX, posOnSurf_strip.y());
292 stripNumber = detEl->
stripNumber(new_position, newId);
294 if (stripNumber < 1) {
296 ATH_MSG_WARNING(
"Position on strip surface = (" << posOnSurf_strip.x() <<
", " << posOnSurf_strip.y() <<
")");
300 newId = idHelper.channelID(layid, idHelper.multilayer(layid),
301 idHelper.gasGap(layid), channelType, stripNumber,
isValid);
302 if(!
isValid && stripNumber != -1) {
308 const double stripHalfPitch = detEl->
channelPitch(newId)*0.5;
329 const double total_charge =
gain*ionized_charge;
334 const double tan_theta = GLODIRE.perp()/GLODIRE.z();
338 const double cluster_posX = posOnSurf_strip.x();
339 double peak_position = CLHEP::RandGaussZiggurat::shoot(condContainers.rndmEngine, cluster_posX,
m_StripResolution*angle_dependency);
347 std::unique_ptr<TF1> clusterProfile = std::make_unique<TF1>(
"fgaus",
348 "[0]*exp(-0.5*(x/[1])^2)+[2]*exp(-0.5*(x/[3])^2)",
356 constexpr
double tolerance_charge = 0.0005;
358 constexpr
unsigned int max_neighbor = 10;
361 for (
unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
362 int currentStrip = stripNumber + iStrip;
363 if (currentStrip > NumberOfStrips)
break;
366 newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid),
367 channelType, currentStrip,
isValid);
376 double x_relative = locpos.x() - peak_position;
378 double charge = clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
380 if (
charge < tolerance_charge)
break;
383 double strip_time = sDigitTimeStrip;
388 int indexFromMiddleStrip = iStrip;
389 if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
390 indexFromMiddleStrip = iStrip - 1;
398 <<
", time = " << strip_time <<
", time offset = " << strip_time-sDigitTimeStrip
399 <<
", neighbor index = " << iStrip
405 for (
unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
406 int currentStrip = stripNumber - iStrip;
407 if (currentStrip < 1)
break;
409 newId = idHelper.channelID(layid,idHelper.multilayer(layid),
410 idHelper.gasGap(layid), channelType, currentStrip,
isValid);
418 double x_relative = locpos.x() - peak_position;
419 double charge = clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
420 if (
charge < tolerance_charge)
break;
423 double strip_time = sDigitTimeStrip;
426 int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
433 <<
", time = " << strip_time <<
", time offset = " << strip_time-sDigitTimeStrip
434 <<
", neighbor index = " << iStrip
449 channelType = sTgcIdHelper::sTgcChannelTypes::Pad;
451 Identifier PAD_ID = idHelper.channelID(layid, idHelper.multilayer(layid),
452 idHelper.gasGap(layid), channelType, 1);
458 const Amg::Vector2D posOnSurf_pad(hitOnSurface_pad.x(), hitOnSurface_pad.y());
464 int padNumber = detEl->
stripNumber(posOnSurf_pad, PAD_ID);
465 if( padNumber == -1 ){
467 const double new_posX = (posOnSurf_pad.x()>0.0)? posOnSurf_pad.x()-length_tolerance
468 : posOnSurf_pad.x()+length_tolerance;
469 const double new_posY = (posOnSurf_pad.y()>0.0)? posOnSurf_pad.y()-length_tolerance
470 : posOnSurf_pad.y()+length_tolerance;
472 padNumber = detEl->
stripNumber(new_position, PAD_ID);
476 ATH_MSG_WARNING(
"Position on pad surface = (" << posOnSurf_pad.x() <<
", " << posOnSurf_pad.y() <<
")");
480 newId = idHelper.channelID(layid, idHelper.multilayer(layid),
481 idHelper.gasGap(layid), channelType, padNumber,
isValid);
499 double deltaX = halfPadWidthX - std::abs(
diff.x());
500 double deltaY= halfPadWidthY - std::abs(
diff.y());
512 int newPhi =
diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
513 int newEta =
diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
517 if (isNeighX && isNeighY && validEta && validPhi){
519 Identifier neigh_ID_X = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
520 channelType, padEtaPhi.first, newPhi,
isValid);
521 Identifier neigh_ID_Y = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
522 channelType, newEta, padEtaPhi.second,
isValid);
523 Identifier neigh_ID_XY = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
524 channelType, newEta, newPhi,
isValid);
529 addDigit(digits, neigh_ID_X, bctag, sDigitTimePad, xQfraction*(1.-yQfraction)*0.5*total_charge);
530 addDigit(digits, neigh_ID_Y, bctag, sDigitTimePad, yQfraction*(1.-xQfraction)*0.5*total_charge);
531 addDigit(digits, neigh_ID_XY, bctag, sDigitTimePad, xQfraction*yQfraction*0.5*total_charge);
532 addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction-yQfraction+xQfraction*yQfraction)*0.5*total_charge);
533 }
else if (isNeighX && validPhi){
535 Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
536 channelType, padEtaPhi.first, newPhi,
isValid);
541 addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction)*0.5*total_charge);
542 addDigit(digits, neigh_ID, bctag, sDigitTimePad, xQfraction*0.5*total_charge);
544 else if (isNeighY && validEta){
546 Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
547 channelType, newEta, padEtaPhi.second,
isValid);
552 addDigit(digits, newId, bctag, sDigitTimePad, (1.-yQfraction)*0.5*total_charge);
553 addDigit(digits, neigh_ID, bctag, sDigitTimePad, yQfraction*0.5*total_charge);
559 addDigit(digits, newId, bctag, sDigitTimePad, 0.5*total_charge);
563 else if(padNumber != -1) {
581 channelType = sTgcIdHelper::sTgcChannelTypes::Wire;
584 Identifier WIREGP_ID = idHelper.channelID(layid, idHelper.multilayer(layid),
585 idHelper.gasGap(layid), channelType, 1);
593 int wiregroupNumber = detEl->
stripNumber(posOnSurf_wire, WIREGP_ID);
594 if( wiregroupNumber == -1 ){
596 const double new_posX = (posOnSurf_wire.x() > 0.0)? posOnSurf_wire.x() - length_tolerance
597 : posOnSurf_wire.x() + length_tolerance;
598 const Amg::Vector2D new_position(new_posX, posOnSurf_wire.y());
599 wiregroupNumber = detEl->
stripNumber(new_position, WIREGP_ID);
601 if (wiregroupNumber < 1) {
603 ATH_MSG_WARNING(
"Position on wire surface = (" << posOnSurf_wire.x() <<
", " << posOnSurf_wire.y() <<
")");
609 newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid),
610 channelType, wiregroupNumber,
isValid);
614 if(wiregroupNumber>=1&&wiregroupNumber<=NumberOfWiregroups) {
615 addDigit(digits, newId, bctag, sDigitTimeWire, total_charge);
618 else if (wiregroupNumber != -1){