ATLAS Offline Software
MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
7 #include "CLHEP/Units/SystemOfUnits.h"
8 #include "CLHEP/Random/RandomEngine.h"
9 #include "CLHEP/Random/RandFlat.h"
10 #include "CLHEP/Random/RandGaussZiggurat.h"
11 #include "CLHEP/Random/RandGamma.h"
13 
14 #include "TF1.h"
15 #include <cmath>
16 #include <iostream>
17 #include <fstream>
18 
19 //---------------------------------------------------
20 // Constructor and Destructor
21 //---------------------------------------------------
22 
25  double meanGasGain,
26  bool doPadChargeSharing)
27  : AthMessaging("sTgcDigitMaker"),
28  m_idHelperSvc{idHelperSvc},
29  m_digitMode(mode),
30  m_meanGasGain{meanGasGain},
31  m_doPadSharing{doPadChargeSharing}{}
32 
34 
35 //------------------------------------------------------
36 // Initialize digitization parameters
37 //------------------------------------------------------
39  // Read arrival time data
41 
42  // Read strip time correction if enabled
43  if (m_doTimeOffsetStrip) {
45  }
46 
47  return StatusCode::SUCCESS;
48 }
49 
50 //------------------------------------------------------
51 // Execute digitization for a given hit
52 //------------------------------------------------------
54  // Extract energy deposit from the hit
55  double energyDeposit = hit->energyDeposit();
56  if (energyDeposit == 0.) return {}; // Ignore hits with no energy deposit
57 
58  // Retrieve the detector element for the given hit
59  const MuonGMR4::MuonDetectorManager* detMgr = condContainers.detMgr;
60 
61  // Convert hit to identifier and retrieve sTgc readout element
62  const Identifier hitId = hit->identify();
63  ATH_MSG_DEBUG("Retrieving detector element for: "<< m_idHelperSvc->toStringDetEl(hitId) << " energyDeposit "<< energyDeposit );
64  const MuonGMR4::sTgcReadoutElement* readoutElement = detMgr->getsTgcReadoutElement(hitId);
65 
66  // Projecting the hit position on wire surface given the hit position
67  // and direction w.r.t. the wire plane
68  Amg::Vector3D locHitDir = xAOD::toEigen(hit->localDirection());
69  Amg::Vector3D locHitPos = xAOD::toEigen(hit->localPosition());
70  ATH_MSG_VERBOSE("sTgc hit: time " << hit->globalTime()
71  << " position " << Amg::toString(locHitPos, 2)
72  << " direction" << Amg::toString(locHitDir, 2)
73  << " mclink " << hit->genParticleLink() << " PDG ID " << hit->pdgId() );
74  ATH_MSG_VERBOSE("Projecting hit to Wire Surface" );
75  const double scale = Amg::intersect<3>(locHitPos, locHitDir, Amg::Vector3D::UnitZ(), 0.).value_or(0);
76  Amg::Vector3D hitOnWireSurf = locHitPos + scale * locHitDir;
77  ATH_MSG_VERBOSE("Hit on Wire Surface: " << Amg::toString(hitOnWireSurf, 2));
78 
79  /* Determine the closest wire and the distance of closest approach
80  * Since most particles pass through the the wire plane between two wires,
81  * the nearest wire should be one of these two wire. Otherwise, the particle's
82  * trajectory is uncommon, and such rare case is not supported yet.
83  *
84  * Finding that nearest wire follows the following steps:
85  * - Compute the distance to the wire at the center of the current wire pitch
86  * - Compute the distance to the other adjacent wire
87  */
88  const sTgcIdHelper& idHelper{m_idHelperSvc->stgcIdHelper()};
89  const Identifier wireLayId = idHelper.channelID(hitId,
90  idHelper.multilayer(hitId),
91  idHelper.gasGap(hitId),
92  ReadoutChannelType::Wire,
93  1);
94 
95  const MuonGMR4::WireGroupDesign& wireDesign{readoutElement->wireDesign(wireLayId)};
96  Amg::Vector2D hitOnWire2D = hitOnWireSurf.block<2,1>(0,0);
97  if(!wireDesign.insideTrapezoid(hitOnWire2D)) {
98  return {};
99  }
100  std::pair<int, int> wireGrpWireNum = wireDesign.wireNumber(hitOnWire2D);
101  int wireNumber = -1;
102  wireNumber = wireDesign.numPitchesToGroup(wireGrpWireNum.first) + wireGrpWireNum.second;
103  // If wire number is invalid, verify if hit is near the edge of the chamber.
104  // Try to get the wire number from the x and y coordinates of the hit position rather than the hit on wire surface.
105  const int numWires = wireDesign.nAllWires();
106  if((wireNumber < 1) || (wireNumber > numWires)) {
107  wireGrpWireNum = wireDesign.wireNumber(locHitPos.block<2,1>(0,0));
108  wireNumber = wireDesign.numPitchesToGroup(wireGrpWireNum.first) + wireGrpWireNum.second;
109  if((wireNumber < 1) || (wireNumber > numWires)) {
110  ATH_MSG_WARNING("Unable to retrieve the wire number, skipping the hit: " << m_idHelperSvc->toString(hitId));
111  return {};
112  }
113  }
114  // Compute the position of the ionization and its distance to the closest wire.
115  Ionization ionization = pointClosestApproach(wireDesign, wireNumber, locHitPos, locHitDir, hit->stepLength());
116  double distToWire = ionization.distance;
117  if(distToWire > 0.) {
118  // Determine on which side of the wire does the particle cross
119  // -1 if particle crosses the wire surface between wireNumber-1 and wireNumber
120  // +1 if particle crosses the wire surface between wireNumber and wireNumber+1
121  int adjacent = 1;
122  if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {adjacent = -1;}
123  // Find the position of the ionization with respect to the adjacent wire
124  Ionization ionizationAdj = pointClosestApproach(wireDesign, wireNumber+adjacent, locHitPos, locHitDir, hit->stepLength());
125  // Determine if the adjacent wire is closer to the hit on wire surface
126  double distToWireAdj = ionizationAdj.distance;
127  if ((distToWireAdj > 0.) && (distToWireAdj < distToWire)) {
128  distToWire = distToWireAdj;
129  wireNumber += adjacent;
130  ionization = std::move(ionizationAdj);
131  }
132  } else {
133  ATH_MSG_DEBUG("Failed to get the distance between the wire number = " << wireNumber
134  << " and hit at " <<Amg::toString(hitOnWireSurf, 2)
135  << ". Number of wires = " << numWires
136  << ", "<<m_idHelperSvc->toStringGasGap(wireLayId));
137  return {};
138  }
139 
140  // Update the position of ionization on the wire surface
141  hitOnWireSurf = ionization.posOnWire;
142  ATH_MSG_VERBOSE("Ionization_info: distance: " << ionization.distance
143  << " posOnTrack: " <<Amg::toString(ionization.posOnSegment, 3)
144  << " posOnWire: " << Amg::toString(ionization.posOnWire, 3)
145  << " hit local Pos: " << Amg::toString(locHitPos,2)
146  << " hit local Dir: " <<Amg::toString(locHitDir, 2)
147  << " EDep: " << hit->energyDeposit() << " EKin: " << hit->kineticEnergy()
148  << " pdgId: " << hit->pdgId()
149  << "gasGap: "<<m_idHelperSvc->toStringGasGap(wireLayId));
150 
151  // Distance should be in the range [0, 0.9] mm, excepting
152  // - particles pass through the wire plane near the edges
153  // - secondary particles created inside the gas gap that go through the gas gap partially.
154  // Most of such particles are not muons and have low kinetic energy.
155  // - particle with trajectory parallel to the sTGC wire plane
156  const double wirePitch = wireDesign.stripPitch();
157  if ((distToWire > 0.) && (std::abs(hit->pdgId()) == 13) && (distToWire > (0.5 * wirePitch))) {
158  ATH_MSG_DEBUG("Distance to the nearest wire (" << distToWire << ") is greater than expected.");
159  ATH_MSG_DEBUG("Hit local Pos: "<<Amg::toString(locHitPos, 2)
160  << " hit local Dir: " <<Amg::toString(locHitDir, 2)
161  << " EDeposited: " << hit->energyDeposit()
162  << " EKinetic: " << hit->kineticEnergy()
163  << " pdgID: " << hit->pdgId()
164  << " gasGap: "<<m_idHelperSvc->toStringGasGap(wireLayId));
165  }
166  // Do not digitize hits that are too far from the nearest wire
167  if (distToWire > wirePitch) {
168  return {};
169  }
170 
171  // Get the gamma pdf parameters associated with the distance of closest approach.
172  const GammaParameter gamParam = getGammaParameter(distToWire);
173  const double par_kappa = gamParam.kParameter;
174  const double par_theta = gamParam.thetaParameter;
175  const double most_prob_time = getMostProbableArrivalTime(distToWire);
176  // Compute the most probable value of the gamma pdf
177  double gamma_mpv = (par_kappa - 1) * par_theta;
178  // If the most probable value is less than zero, then set it to zero
179  if (gamma_mpv < 0.) {gamma_mpv = 0.;}
180  const double t0_par = most_prob_time - gamma_mpv;
181 
182  // Digit time follows a gamma distribution, so a value val is
183  // chosen using a gamma random generator then is shifted by t0
184  // to account for drift time.
185  // Note: CLHEP::RandGamma takes the parameters k and lambda,
186  // where lambda = 1 / theta.
187  double digitTime = t0_par + CLHEP::RandGamma::shoot(condContainers.rndEngine, par_kappa, 1/par_theta);
188 
189  // Sometimes, digitTime is negative because t0_par can be negative.
190  // In such case, discard the negative value and shoot RandGamma for another value.
191  // However, if that has already been done many times then set digitTime to zero
192  // in order to avoid runaway loop.
193 
194  constexpr int shoot_limit = 4;
195  int shoot_counter = 0;
196  while (digitTime < 0.) {
197  if (shoot_counter > shoot_limit) {
198  digitTime = 0.;
199  break;
200  }
201  digitTime = t0_par + CLHEP::RandGamma::shoot(condContainers.rndEngine, par_kappa, 1/par_theta);
202  ++shoot_counter;
203  }
204 
205  ATH_MSG_DEBUG("sTgcDigitMaker distance = " << distToWire
206  << ", time = " << digitTime
207  << ", k parameter = " << par_kappa
208  << ", theta parameter = " << par_theta
209  << ", most probable time = " << most_prob_time);
210 
212  if (condContainers.efficiencies) {
213  const double efficiency = condContainers.efficiencies->getEfficiency(wireLayId);
214  // Lose Hits to match HV efficiency
215  if (CLHEP::RandFlat::shoot(condContainers.rndEngine,0.0,1.0) > efficiency) return {};
216  }
217 
218  sTgcDigitVec digits{};
219 
220  double sDigitTimeWire = digitTime;
221  double sDigitTimePad = sDigitTimeWire;
222  double sDigitTimeStrip = sDigitTimeWire;
223 
224  uint16_t bctag = 0;
225 
226  //##################################################################################
227  //######################################### strip readout ##########################
228  //##################################################################################
229  ATH_MSG_DEBUG("sTgcDigitMaker::strip response ");
230  bool isValid = false;
231  int channelType = ReadoutChannelType::Strip;
232  const Identifier stripLayId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId),
233  idHelper.gasGap(wireLayId), channelType, 1, isValid);
234 
235  const MuonGMR4::StripDesign& stripDesign{readoutElement->stripDesign(stripLayId)};
236  // Rotate clockwise 90 in xy plane to get to strip plane orientation, and then extract x,y
237  const Amg::Vector2D hitOnStripSurf = (Amg::getRotateZ3D(-90. * Gaudi::Units::deg) * hitOnWireSurf).block<2,1>(0,0);
238 
239  bool insideBounds = stripDesign.insideTrapezoid(hitOnStripSurf);
240  if(!insideBounds) {
241  ATH_MSG_DEBUG("Outside of the strip surface boundary : " << m_idHelperSvc->toString(stripLayId) << "; local position " <<Amg::toString(hitOnStripSurf, 2));
242  return {};
243  }
244 
245  //************************************ find the nearest readout element **************************************
246  // Required precision on length in mm
247  constexpr double tolerance_length = 0.01*Gaudi::Units::mm;
248  int stripNumber = stripDesign.stripNumber(hitOnStripSurf);
249  if( stripNumber < 0){
250  // Verify if the energy deposit is at the boundary
251  const double newPosX = (hitOnStripSurf.x() > 0.0)? hitOnStripSurf.x() - tolerance_length
252  : hitOnStripSurf.x() + tolerance_length;
253  const Amg::Vector2D newPos(newPosX, hitOnStripSurf.y());
254  stripNumber = stripDesign.stripNumber(newPos);
255  // Skip hit if still unable to obtain strip number
256  if (stripNumber < 0) {
257  ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperSvc->toString(stripLayId) );
258  ATH_MSG_WARNING("Position on strip surface = (" << hitOnStripSurf.x() << ", " << hitOnStripSurf.y() << ")");
259  return {};
260  }
261  }
262  isValid = false;
263  const Identifier stripHitId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId),
264  idHelper.gasGap(stripLayId), channelType, stripNumber, isValid);
265  if(!isValid && stripNumber != -1) {
266  ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(stripHitId) );
267  return {};
268  }
269 
270  const int numStrips = stripDesign.numStrips();
271  const double stripHalfPitch = 0.5 * stripDesign.stripPitch();
272 
273  //************************************ conversion of energy to charge **************************************
274 
275  // Typical ionized charge in pC per keV deposited. The constant is determined from ionization
276  // study with Garfield program. A note titled "Charge Energy Relation" which outlines
277  // conversion can be found here:
278  // https://gitlab.cern.ch/carleton-ATLAS/carleton-nsw/stgc-ionization-gain/-/blob/master/Charge_Energy_Relation.pdf?ref_type=heads
279  const double ionized_charge = (5.65E-6)*energyDeposit/CLHEP::keV;
280 
281  // To get avalanche gain, polya function is taken from Blum paper https://inspirehep.net/literature/807304
282  // m_polyaFunction = new TF1("m_polyaFunction","(1.0/[1])*(TMath::Power([0]+1,[0]+1)/TMath::Gamma([0]+1))*TMath::Power(x/[1],[0])*TMath::Exp(-([0]+1)*x/[1])",0,3000000);
283 
284  // Mean value for total gain due to E field;
285  // To calculate this gain from polya distibution, we replace in gamma PDF:
286  // alpha = 1+theta and
287  // beta = 1+theta/mean
288  // With these substitutions, gamma PDF gives the same sampling values as those from polya PDF.
289  const double gain = CLHEP::RandGamma::shoot(condContainers.rndEngine, 1. + m_theta, (1. + m_theta)/m_meanGasGain);
290 
291  // total charge after avalanche
292  const double total_charge = gain*ionized_charge;
293 
294  // Charge Spread including tan(theta) resolution term.
295  const double tan_theta = locHitDir.perp()/locHitDir.z();
296  // The angle dependance on strip resolution goes as tan^2(angle)
297  const double angle_dependency = std::hypot(m_posResIncident, m_posResAngular * tan_theta);
298 
299  const double cluster_posX = hitOnStripSurf.x();
300  double peak_position = CLHEP::RandGaussZiggurat::shoot(condContainers.rndEngine, cluster_posX, m_StripResolution*angle_dependency);
301 
302  // Each readout plane reads about half the total charge produced on the wire,
303  // including a tan(theta) term to describe the increase of charge with incident angle
304  const double norm = 0.5 * total_charge;
305 
306  // Lower limit on strip charge (arbitrary limit), in pC, which has the same units as the parameter ionized_charge.
307  constexpr double tolerance_charge = 0.0005;
308  // Set a maximum number of neighbour strips to avoid very long loop. This is an arbitrary number.
309  constexpr unsigned int max_neighbor = 10;
310 
311  // Spread charge on the strips that are on the upper half of the strip cluster
312  for (unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
313  int currentStrip = stripNumber + iStrip;
314  if (currentStrip > numStrips) {
315  ATH_MSG_VERBOSE("Breaking the upper half strip loop stripNumber: " << currentStrip << " > " << numStrips);
316  break;
317  }
318 
319  // Get the strip identifier and create the digit
320  isValid = false;
321  const Identifier currentStripId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId), idHelper.gasGap(stripLayId),
322  channelType, currentStrip, isValid);
323  if (isValid) {
324  Amg::Vector2D currentStripPos = readoutElement->localChannelPosition(currentStripId);
325  if (currentStripPos == Amg::Vector2D::Zero()) {
326  ATH_MSG_VERBOSE("Failed to obtain local position for identifier " << m_idHelperSvc->toString(currentStripId) );
327  }
328 
329  // Estimate the digit charge
330  // Position with respect to the peak of the charge curve
331  double x_relative = currentStripPos.x() - peak_position;
332  // In clusterProfile curve, position should be in the units of strip channel
333  double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * norm * chargeIntegral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
334  // If charge is too small, stop creating neighbor strip
335  if (charge < tolerance_charge) {
336  ATH_MSG_VERBOSE("Breaking the upper half strip loop stripCharge: " << charge << " < " << tolerance_charge);
337  break;
338  }
339 
340  // Estimate digit time
341  double strip_time = sDigitTimeStrip;
342  // Strip time response can be delayed due to the resistive layer.
343  // A correction would be required if the actual VMM front-end doesn't re-align the strip timing.
344  if (m_doTimeOffsetStrip) {
345  // Determine how far the current strip is from the middle strip
346  int indexFromMiddleStrip = iStrip;
347  if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
348  indexFromMiddleStrip = iStrip - 1;
349  // Add time delay due to resistive layer
350  strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
351  }
352 
353  addDigit(digits, currentStripId, bctag, strip_time, charge);
354 
355  ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
356  << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
357  << ", neighbor index = " << iStrip
358  << ", strip position = "<<Amg::toString(currentStripPos, 2));
359  }
360  }
361 
362  // The lower half of the strip cluster
363  for (unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
364  int currentStrip = stripNumber - iStrip;
365  if (currentStrip < 1) {
366  ATH_MSG_VERBOSE("Breaking the lower half strip loop stripNumber: " << currentStrip << " < 1 ");
367  break;
368  }
369 
370  isValid = false;
371  const Identifier currentStripId = idHelper.channelID(stripLayId, idHelper.multilayer(stripLayId), idHelper.gasGap(stripLayId),
372  channelType, currentStrip, isValid);
373  if (isValid) {
374  Amg::Vector2D currentStripPos = readoutElement->localChannelPosition(currentStripId);
375  if (currentStripPos == Amg::Vector2D::Zero()) {
376  ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(currentStripId) );
377  }
378 
379  // Estimate the digit charge
380  double x_relative = currentStripPos.x() - peak_position;
381  double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * norm * chargeIntegral((x_relative/(2*stripHalfPitch) - 0.5), (x_relative/(2*stripHalfPitch) + 0.5));
382  if (charge < tolerance_charge) {
383  ATH_MSG_VERBOSE("Breaking the lower half strip loop stripCharge: " << charge << " < " << tolerance_charge);
384  break;
385  }
386 
387  // Estimate digit time
388  double strip_time = sDigitTimeStrip;
389  // Time delay due to resistive layer
390  if (m_doTimeOffsetStrip) {
391  int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
392  strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
393  }
394 
395  addDigit(digits, currentStripId, bctag, strip_time, charge);
396 
397  ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
398  << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
399  << ", neighbor index = " << iStrip
400  << ", strip position = " << Amg::toString(currentStripPos, 2));
401  }
402  }
403  // end of strip digitization
404 
405  if(m_digitMode == digitMode::StripsOnly) {
406  ATH_MSG_WARNING("Only digitize strip response !");
407  return digits;
408  }
409 
410  //##################################################################################
411  //######################################### pad readout ##########################
412  //##################################################################################
413  ATH_MSG_DEBUG("sTgcDigitMaker::pad response ");
414  channelType = ReadoutChannelType::Pad;
415 
416  Identifier padLayId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId),
417  idHelper.gasGap(wireLayId), channelType, 1);
418 
419  const MuonGMR4::PadDesign& padDesign{readoutElement->padDesign(padLayId)};
420  //************************************ find the nearest readout element **************************************
421 
422  // Since the orientation of pad layer is the same as wirelayer, the x,y can be taken directly
423  // from the hit position on wire surface
424  const Amg::Vector2D hitOnPadSurf = hitOnWireSurf.block<2,1>(0,0);
425 
426 
427  insideBounds = padDesign.insideTrapezoid(hitOnPadSurf);
428 
429  if(insideBounds) {
430  std::pair<uint, uint> padEtaPhi = padDesign.channelNumber(hitOnPadSurf);
431  unsigned int padEta = padEtaPhi.first;
432  unsigned int padPhi = padEtaPhi.second;
433  if( padEta < 1 || padPhi < 1){
434  // Verify if the energy deposit is at the boundary
435  const double newPosX = (hitOnPadSurf.x()>0.0)? hitOnPadSurf.x() - tolerance_length
436  : hitOnPadSurf.x() + tolerance_length;
437  const double newPosY = (hitOnPadSurf.y()>0.0)? hitOnPadSurf.y() - tolerance_length
438  : hitOnPadSurf.y() + tolerance_length;
439  const Amg::Vector2D newPosition(newPosX, newPosY);
440  padEtaPhi = padDesign.channelNumber(newPosition);
441  padEta = padEtaPhi.first;
442  padPhi = padEtaPhi.second;
443  // Skip hit if still unable to obtain pad number
444  if( padEta < 1 || padPhi < 1){
445  ATH_MSG_WARNING("Failed to obtain pad number for " << m_idHelperSvc->toString(padLayId) );
446  ATH_MSG_WARNING("Position on pad surface = (" << hitOnPadSurf.x() << ", " << hitOnPadSurf.y() << ")");
447  return digits;
448  }
449  }
450  isValid = false;
451  const Identifier padHitId = idHelper.padID(padLayId, idHelper.multilayer(padLayId),
452  idHelper.gasGap(padLayId), channelType, padEta, padPhi, isValid);
453  if(isValid) {
454  // Find centre position of pad
455  Amg::Vector2D padPos = readoutElement->localChannelPosition(padHitId);
456  if (padPos == Amg::Vector2D::Zero()) {
457  ATH_MSG_ERROR("Failed to obtain local position for neighbor pad with identifier " << m_idHelperSvc->toString(padHitId) );
458  return digits;
459  }
460 
461  // Pad sharing needs to look at position on hit vs pad boundaries
462  const Amg::Vector2D diff = hitOnPadSurf - padPos;
463  double halfPadHeight = 0.5 * readoutElement->padHeight(padHitId);
464  // Extracting pad corners to get the width of the pad at hit position
465  std::array<Amg::Vector2D, 4> padHitCorners = padDesign.padCorners(padEtaPhi);
466  double padBottomBase = (padHitCorners[0] - padHitCorners[1]).norm();
467  double padTopBase = (padHitCorners[2] - padHitCorners[3]).norm();
468  // Linearly interpolating pad width (as done by channelWidth function in R3)
469  // Interpolating pad width at the y-position of the hit (diff.y()) using pad top and bottom bases.
470  // The expression below computes the width at a relative vertical position inside the pad.
471  // Specifically: (1 + diff.y() / halfPadHeight) * 0.25 = (y_from_bottom / padHeight)
472  double halfPadWidth = 0.5 * padBottomBase + 0.25 * (1 + diff.y()/halfPadHeight) * (padTopBase - padBottomBase);
473 
474  // Charge sharing happens within 4mm window for pads
475  // i.e. the charge is spread over a 2D gaussian of total radius about 4mm
476  // This value depends on actual width of the distribution, but we don't have the
477  // width of the charge distribution for pads. So lets consider 2.5*sigma, where
478  // sigma is the width of the charge distribution for strips.
479  double deltaX = halfPadWidth - std::abs(diff.x());
480  double deltaY = halfPadHeight - std::abs(diff.y());
481  bool isNeighX = deltaX < 2.5*m_clusterParams[1];
482  bool isNeighY = deltaY < 2.5*m_clusterParams[1];
483  // Pad width can be calculated to be very slightly larger than it should due to rounding errors
484  // So if a hit falls on a given pad but is "outside" the width, just define it to be on the boundary of 2 pads.
485  if (deltaX < 0.) deltaX = 0.1;
486  if (deltaY < 0.) deltaY = 0.1;
487 
488  if (m_doPadSharing && (isNeighX || isNeighY)){
489  // Phi == 1 at the right in the local geometry
490  // In local coordinates, the pad to the right has phi' = phi - 1
491  unsigned int newPhi = diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
492  unsigned int newEta = diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
493  bool validEta = newEta > 0 && newEta < readoutElement->numPadEta(padHitId) + 1;
494  bool validPhi = newPhi > 0 && newPhi < readoutElement->numPadPhi(padHitId) + 1;
495 
496  if (isNeighX && isNeighY && validEta && validPhi){
497  // 4 pads total; makes life a bit harder. Corner of 4 valid pads
498  isValid = false;
499  const Identifier neigh_ID_X = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
500  channelType, padEta, newPhi, isValid);
501  isValid = false;
502  const Identifier neigh_ID_Y = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
503  channelType, newEta, padPhi, isValid);
504  isValid = false;
505  const Identifier neigh_ID_XY = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
506  channelType, newEta, newPhi, isValid);
507  double xQfraction = getPadChargeFraction(deltaX);
508  double yQfraction = getPadChargeFraction(deltaY);
509 
510  // Main pad gets 1 - Qfraction of total
511  addDigit(digits, neigh_ID_X, bctag, sDigitTimePad, xQfraction*(1.-yQfraction)*0.5*total_charge);
512  addDigit(digits, neigh_ID_Y, bctag, sDigitTimePad, yQfraction*(1.-xQfraction)*0.5*total_charge);
513  addDigit(digits, neigh_ID_XY, bctag, sDigitTimePad, xQfraction*yQfraction*0.5*total_charge);
514  addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-xQfraction-yQfraction+xQfraction*yQfraction)*0.5*total_charge);
515  } else if (isNeighX && validPhi){
516  // There is only 1 neighbor, immediately to the left or right.
517  isValid = false;
518  const Identifier neigh_ID = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
519  channelType, padEta, newPhi, isValid);
520 
521  // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
522  // Main pad gets 1 - Qfraction of total
523  double xQfraction = getPadChargeFraction(deltaX);
524  addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-xQfraction)*0.5*total_charge);
525  addDigit(digits, neigh_ID, bctag, sDigitTimePad, xQfraction*0.5*total_charge);
526  }
527  else if (isNeighY && validEta){
528  // There is only 1 neighbor, immediately above or below
529  isValid = false;
530  const Identifier neigh_ID = idHelper.padID(padHitId, idHelper.multilayer(padHitId), idHelper.gasGap(padHitId),
531  channelType, newEta, padPhi, isValid);
532 
533  // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
534  // Main pad gets 1 - Qfraction of total
535  double yQfraction = getPadChargeFraction(deltaX);
536  addDigit(digits, padHitId, bctag, sDigitTimePad, (1.-yQfraction)*0.5*total_charge);
537  addDigit(digits, neigh_ID, bctag, sDigitTimePad, yQfraction*0.5*total_charge);
538  }
539 
540  }
541  else{
542  // No charge sharing: hit is nicely isolated within pad
543  addDigit(digits, padHitId, bctag, sDigitTimePad, 0.5*total_charge);
544  }
545 
546  }
547  else if(padEta != 0 || padPhi != 0) {
548  ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(padHitId) );
549  }
550  }
551  else {
552  ATH_MSG_DEBUG("Outside of the pad surface boundary :" << m_idHelperSvc->toString(padLayId)<< " local position " <<Amg::toString(hitOnPadSurf, 2));
553  }
554 
555  if(m_digitMode == digitMode::StripsAndPads) {
556  ATH_MSG_WARNING("Only digitize strip/pad response !");
557  return digits;
558  }
559 
560  //##################################################################################
561  //######################################### wire readout ##########################
562  //##################################################################################
563  ATH_MSG_DEBUG("sTgcDigitMaker::wire response ");
564  channelType = ReadoutChannelType::Wire;
565 
566  // wireLayId already defined at the start
567  //************************************ find the nearest readout element **************************************
568  insideBounds = wireDesign.insideTrapezoid(hitOnWireSurf.block<2,1>(0,0));
569 
570  if(insideBounds) {
571  // Determine the wire number
572  int wiregroupNumber = (wireDesign.wireNumber(hitOnWireSurf.block<2,1>(0,0))).first;
573  if( wiregroupNumber < 1){
574  // Verify if the energy deposit is at the boundary
575  const double newPosX = (hitOnWireSurf.x() > 0.0)? hitOnWireSurf.x() - tolerance_length
576  : hitOnWireSurf.x() + tolerance_length;
577  const Amg::Vector2D newPos(newPosX, hitOnWireSurf.y());
578  wiregroupNumber = (wireDesign.wireNumber(newPos)).first;
579  // Skip hit if still unable to obtain wire number
580  if (wiregroupNumber < 1) {
581  ATH_MSG_WARNING("Failed to obtain wire number " << m_idHelperSvc->toString(wireLayId) );
582  ATH_MSG_WARNING("Position on wire surface = (" << hitOnWireSurf.x() << ", " << hitOnWireSurf.y() << ")");
583  return digits;
584  }
585  }
586 
587  // Find ID of the actual wiregroup
588  isValid = false;
589  const Identifier wireGroupHitId = idHelper.channelID(wireLayId, idHelper.multilayer(wireLayId), idHelper.gasGap(wireLayId),
590  channelType, wiregroupNumber, isValid);
591 
592  if(isValid) {
593  int numWireGroups = readoutElement->numChannels(wireLayId);
594  if(wiregroupNumber >=1 && wiregroupNumber <= numWireGroups) {
595  addDigit(digits, wireGroupHitId, bctag, sDigitTimeWire, total_charge);
596  }
597  } // end of if(isValid)
598  else if (wiregroupNumber != -1){
599  ATH_MSG_ERROR("Failed to obtain wiregroup identifier " << m_idHelperSvc->toString(wireGroupHitId));
600  }
601  }
602  else {
603  ATH_MSG_DEBUG("Outside of the wire surface boundary :" << m_idHelperSvc->toString(wireLayId)<< " local position " <<Amg::toString(hitOnWireSurf, 2));
604  }
605  // end of wire digitization
606 
607  return digits;
608 }
609 
610 //------------------------------------------------------
611 // Compute closest approach between hit trajectory and wire
612 //------------------------------------------------------
614  int wireNumber,
615  const Amg::Vector3D& locHitPos,
616  const Amg::Vector3D& locHitDir,
617  const double stepLength) const {
618 
619  constexpr double angular_tolerance = 1e-3;
620  // Position of the ionization
621  Ionization ionization;
622  // Finding smallest distance and the points at the smallest distance.
623  // The smallest distance between two lines is perpendicular to both lines.
624  // Previous logic was to find the perpendicular distance between two lines using projection geometry
625  // We can construct two lines in the wire surface local coordinate frame:
626  // - one for the hit segment with equation h0 + t * v_h, where h0 is a point
627  // and v_h is the unit vector of the hit segment
628  // - another for the wire with similar equation w0 + s * v_w, where w0 is a
629  // point and v_w is the unit vector of the wire line
630  // Then it is possible to determine the closest points on each line
631  // by requiring that the vector between them is perpendicular to both:
632  // 1. (h0 + t*v_h - w0 - s*v_w) · v_h = 0
633  // 2. (h0 + t*v_h - w0 - s*v_w) · v_w = 0
634 
635  // We have replaced this logic with using Amg::intersect<3> method
636 
637  // Geometry setup
638  const double wirePitch = wireDesign.stripPitch();
639  const double wirePosX = wireDesign.firstStripPos().x() + (wireNumber - 1) * wirePitch;
640  const Amg::Vector3D wireDir(wireDesign.stripDir().x(), wireDesign.stripDir().y(), 0.0);
641  const Amg::Vector3D wirePos(wirePosX, locHitPos.y(), 0.);
642 
643  // Use Amg::intersect to find closest point on hit segment to wire plane
644  std::optional<double> scaleHit = Amg::intersect<3>(locHitPos, locHitDir, Amg::Vector3D::UnitZ(), 0);
645  if (!scaleHit || std::abs(std::abs(wireDir.dot(locHitDir)) - 1.0) < angular_tolerance) {
646  ATH_MSG_DEBUG("The track segment is parallel to the wire, position of digit is undefined");
647  ionization.posOnSegment = locHitPos;
648  ionization.posOnWire = wirePos;
649  ionization.distance = std::hypot(locHitPos.x() - wirePosX, locHitPos.z());
650  return ionization;
651  }
652  // Position on hit segment
653  Amg::Vector3D ionizationPos = locHitPos + scaleHit.value() * locHitDir;
654 
655  if (scaleHit.value() > stepLength) {
656  ionization.posOnSegment = locHitPos;
657  const Amg::Vector3D closestPointToWirePlane = locHitPos + stepLength * locHitDir;
658  ionization.posOnWire = wirePos;
659  ionization.distance = std::abs(closestPointToWirePlane.z());
660  return ionization;
661  }
662 
663  // Project ionization position onto the wire line
664  const double scaleWire = (ionizationPos - wirePos).dot(wireDir);
665  const Amg::Vector3D closestPointOnWire = wirePos + scaleWire * wireDir;
666 
667  // Fill ionization result
668  ionization.posOnSegment = ionizationPos;
669  ionization.posOnWire = closestPointOnWire;
670  ionization.distance = (ionizationPos - closestPointOnWire).mag();
671 
672  return ionization;
673 }
674 
675 //------------------------------------------------------
676 // Adds a digit to the appropriate cache
677 //------------------------------------------------------
678 void sTgcDigitMaker::addDigit(sTgcDigitVec& digits,
679  const Identifier& id,
680  const uint16_t bctag,
681  const double digittime,
682  const double charge) {
683 
684  constexpr double tolerance = 0.1;
685  if (!std::ranges::any_of(digits, [&](std::unique_ptr<sTgcDigit>& known) {
686  return known->identify() == id && std::abs(digittime - known->time()) < tolerance;
687  })) {
688  digits.push_back(std::make_unique<sTgcDigit>(id, bctag, digittime, charge, 0, 0));
689  }
690 }
691 
692 //------------------------------------------------------
693 // Reads time arrival data file
694 //------------------------------------------------------
696  const std::string file_name = "sTGC_Digitization_timeArrival.dat";
697  std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
698  if(file_path.empty()) {
699  ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
700  return StatusCode::FAILURE;
701  }
702 
703  std::ifstream ifs{file_path, std::ios::in};
704  if(ifs.bad()) {
705  ATH_MSG_FATAL("sTgcDigitMaker: Failed to open time of arrival file " << file_name );
706  return StatusCode::FAILURE;
707  }
708 
709  // Read the sTGC_Digitization_timeWindowOffset.dat file
710  std::string line;
711  GammaParameter param{};
712  while (std::getline(ifs, line)) {
713  std::string key;
714  std::istringstream iss(line);
715  iss >> key;
716  if (key.compare("bin") == 0) {
717  iss >> param.lowEdge >> param.kParameter >> param.thetaParameter;
718  m_gammaParameter.push_back(param);
719  } else if (key.compare("mpv") == 0) {
720  double mpt;
721  int idx{0};
722  while (iss >> mpt) {m_mostProbableArrivalTime[idx++] = mpt;}
723  }
724  }
725  ifs.close();
726  return StatusCode::SUCCESS;
727 }
728 
729 //------------------------------------------------------
730 // Retrieves gamma distribution parameters based on distance
731 //------------------------------------------------------
733  const double d = std::abs(distance);
734  // Find the parameters assuming the container is sorted in ascending order of 'lowEdge'
735  if (d < m_gammaParameter.front().lowEdge) {
736  return m_gammaParameter.front();
737  }
738  int index{-1};
739  for (const auto& par: m_gammaParameter) {
740  if (d < par.lowEdge) {
741  break;
742  }
743  ++index;
744  }
745  return m_gammaParameter.at(index);
746 }
747 
748 //------------------------------------------------------
749 // Computes the most probable arrival time based on the distance of closest approach
750 //------------------------------------------------------
752  const double d{std::abs(distance)};
753  double mpt{0};
754  for (size_t t = 0 ; t < m_mostProbableArrivalTime.size(); ++t){
755  mpt += m_mostProbableArrivalTime[t] * std::pow(d, t);
756  }
757  return mpt;
758 }
759 //------------------------------------------------------
760 // Computes the charge on each strip
761 //------------------------------------------------------
762 
763 double sTgcDigitMaker::chargeIntegral(double N, double M) const {
764 
765  double term1 = 0.25 * std::erf( M / (std::sqrt(2) * m_clusterParams[0]));
766  double term2 = 0.25 * std::erf( N / (std::sqrt(2) * m_clusterParams[0]));
767  double term3 = 0.25 * std::erf( M / (std::sqrt(2) * m_clusterParams[1]));
768  double term4 = 0.25 * std::erf( N / (std::sqrt(2) * m_clusterParams[1]));
769 
770  return (term1 - term2 + term3 - term4);
771 }
772 //------------------------------------------------------
773 // Reads strip time offset data file
774 //------------------------------------------------------
776  const std::string file_name = "sTGC_Digitization_timeOffsetStrip.dat";
777  std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
778  if(file_path.empty()) {
779  ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
780  return StatusCode::FAILURE;
781  }
782 
783  // Open the sTGC_Digitization_timeOffsetStrip.dat file
784  std::ifstream ifs{file_path, std::ios::in};
785  if(ifs.bad()) {
786  ATH_MSG_FATAL("Failed to open time of arrival file " << file_name );
787  return StatusCode::FAILURE;
788  }
789 
790  // Initialize the container to store the time offset.
791  // The number of parameters, 6, corresponds to the number of lines to be read
792  // from sTGC_Digitization_timeOffsetStrip.dat.
793  // Setting the default offset to 0 ns.
794  std::string line;
795  size_t index{0};
796  double value{0.0};
797  while (std::getline(ifs, line)) {
798  std::string key;
799  std::istringstream iss(line);
800  iss >> key;
801  if (key.compare("strip") == 0) {
802  iss >> index >> value;
803  if (index >= m_timeOffsetStrip.size()) continue;
805  }
806  }
807  return StatusCode::SUCCESS;
808 }
809 
810 //------------------------------------------------------
811 // Retrieves time offset for strip clusters
812 //------------------------------------------------------
813 double sTgcDigitMaker::getTimeOffsetStrip(size_t neighbor_index) const {
814  return m_timeOffsetStrip.at(std::min(neighbor_index, m_timeOffsetStrip.size() -1));
815 }
816 
817 //------------------------------------------------------
818 // Computes pad charge fraction
819 //------------------------------------------------------
821  // The charge fraction that is found past a distance x away from the
822  // centre of a 2D gaussian distribution of width of cluster profile is
823  // described by a modified error function.
824 
825  // The modified error function perfectly describes
826  // the pad charge sharing distribution figure 16 of the sTGC
827  // testbeam paper https://arxiv.org/pdf/1509.06329.pdf
828  return 0.5 * (1.0 - std::erf( distance / (std::sqrt(2) * m_clusterParams[1])));
829 }
MuonGMR4::sTgcReadoutElement::localChannelPosition
Amg::Vector2D localChannelPosition(const Identifier &measId) const
Returns the local pad/strip/wireGroup position.
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
MuonGMR4::sTgcReadoutElement::padHeight
double padHeight(const Identifier &measId) const
Returns the height of all the pads that are not adjacent to the bottom edge of the trapezoid active a...
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
MuonGMR4::StripDesign
Definition: StripDesign.h:30
Muon::IMuonIdHelperSvc::stgcIdHelper
virtual const sTgcIdHelper & stgcIdHelper() const =0
access to TgcIdHelper
sTgcDigitMaker::m_posResIncident
double m_posResIncident
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:158
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
MuonGMR4::WireGroupDesign
Definition: WireGroupDesign.h:23
AthCheckMacros.h
MuonGMR4::MuonDetectorManager
Definition: MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/MuonDetectorManager.h:62
MuonGMR4::sTgcReadoutElement::numPadEta
unsigned int numPadEta(const Identifier &measId) const
Returns the number of pads in the eta direction in the given layer.
MuonGMR4::sTgcReadoutElement::numChannels
unsigned int numChannels(const Identifier &measId) const
Returns the number of strips / wires / pads in a given gasGap.
makeComparison.deltaY
int deltaY
Definition: makeComparison.py:44
sTgcDigitMaker::getTimeOffsetStrip
double getTimeOffsetStrip(size_t neighbor_index) const
Get digit time offset of a strip depending on its relative position to the strip at the centre of the...
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:860
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
index
Definition: index.py:1
sTgcDigitMaker::readFileOfTimeOffsetStrip
StatusCode readFileOfTimeOffsetStrip()
Read share/sTGC_Digitization_timeOffsetStrip.dat.
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:821
sTgcDigitMaker::m_timeOffsetStrip
std::array< double, 6 > m_timeOffsetStrip
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:137
hist_file_dump.d
d
Definition: hist_file_dump.py:142
sTgcDigitMaker::m_clusterParams
static constexpr std::array< double, 2 > m_clusterParams
Definition: MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h:175
sTgcDigitMaker::Ionization
Ionization object with distance, position on the hit segment and position on the wire.
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:84
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
sTgcDigitMaker::chargeIntegral
double chargeIntegral(double N, double M) const
Definition: MuonPhaseII/MuonDigitization/sTgcDigitizationR4/src/sTgcDigitMaker.cxx:763
MuonGMR4::sTgcReadoutElement::stripDesign
const StripDesign & stripDesign(const Identifier &measId) const
Retrieves the readoutElement Layer given the Identifier/Hash.
MuonGMR4::sTgcReadoutElement::numPadPhi
unsigned int numPadPhi(const Identifier &measId) const
Returns the number of pads in the Phi direction in the given gasGap layer.
CaloCondBlobAlgs_fillNoiseFromASCII.gain
gain
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:109
sTgcDigitMaker::getMostProbableArrivalTime
double getMostProbableArrivalTime(double distance) const
Get the most probable time of arrival.
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:810
sTgcDigitMaker::m_chargeAngularFactor
double m_chargeAngularFactor
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:164
InDet::adjacent
bool adjacent(unsigned int strip1, unsigned int strip2)
Definition: SCT_ClusteringTool.cxx:45
deg
#define deg
Definition: SbPolyhedron.cxx:17
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
sTgcDigitMaker::executeDigi
sTgcDigitVec executeDigi(const DigiConditions &condContainers, const sTGCSimHit &hit) const
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:72
athena.value
value
Definition: athena.py:124
Muon::IMuonIdHelperSvc::toStringDetEl
virtual std::string toStringDetEl(const Identifier &id) const =0
print all fields up to detector element to string
sTgcDigitMaker::readFileOfTimeArrival
StatusCode readFileOfTimeArrival()
Read share/sTGC_Digitization_timeArrival.dat.
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:755
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
sTgcDigitMaker::DigiConditions
Digitize a given hit, determining the time and charge spread on wires, pads and strips.
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:58
TimedHitPtr< xAOD::MuonSimHit >
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
sTgcDigitMaker::GammaParameter
Parameters of a gamma probability distribution function, required for estimating wire digit's time of...
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:75
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
isValid
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
Definition: AtlasPID.h:867
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
sTgcDigitMaker::~sTgcDigitMaker
virtual ~sTgcDigitMaker()
dq_defect_bulk_create_defects.line
line
Definition: dq_defect_bulk_create_defects.py:27
sTgcDigitMaker::GammaParameter::thetaParameter
double thetaParameter
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:78
Trk::energyDeposit
@ energyDeposit
Definition: MeasurementType.h:32
Amg::getRotateZ3D
Amg::Transform3D getRotateZ3D(double angle)
get a rotation transformation around Z-axis
Definition: GeoPrimitivesHelpers.h:270
sTgcDigitMaker::pointClosestApproach
Ionization pointClosestApproach(const MuonGM::sTgcReadoutElement *readoutEle, const Identifier &id, int wireNumber, const Amg::Vector3D &preStepPos, const Amg::Vector3D &postStepPos) const
Determine the points where the distance between two segments is smallest.
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:642
physics_parameters.file_name
string file_name
Definition: physics_parameters.py:32
MuonGMR4::StripDesign::stripDir
const Amg::Vector2D & stripDir() const
Vector pointing along the strip.
sTgcDigitMaker::m_digitMode
digitMode m_digitMode
define offsets and widths of time windows for signals from wiregroups and strips.
Definition: MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h:185
MuonGMR4::sTgcReadoutElement::wireDesign
const WireGroupDesign & wireDesign(const Identifier &measId) const
Retrieves the readoutElement Layer given the Identifier/Hash.
efficiency
void efficiency(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:128
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
sTgcDigitMaker::m_posResAngular
double m_posResAngular
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:159
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
xAOD::uint16_t
setWord1 uint16_t
Definition: eFexEMRoI_v1.cxx:93
sTgcDigitMaker::getGammaParameter
GammaParameter getGammaParameter(double distance) const
Find the gamma pdf parameters of a given distance.
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:795
MuonGMR4::WireGroupDesign::wireNumber
std::pair< int, int > wireNumber(const Amg::Vector2D &extPos) const
Returns a pair where the first component indicate the wire group number and the second one returns th...
MuonGMR4::StripDesign::insideTrapezoid
bool insideTrapezoid(const Amg::Vector2D &extPos) const
Checks whether an external point is inside the trapezoidal area.
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Muon::DigitEffiData::getEfficiency
double getEfficiency(const Identifier &channelId, bool isInnerQ1=false) const
Returns the signal generation efficiency of the sTgc channel.
Definition: DigitEffiData.cxx:38
Muon::IMuonIdHelperSvc::toStringGasGap
virtual std::string toStringGasGap(const Identifier &id) const =0
print all fields up to gas gap to string
athena.file_path
file_path
Definition: athena.py:94
sTgcDigitMaker::m_gammaParameter
std::vector< GammaParameter > m_gammaParameter
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:132
AthMessaging
Class to provide easy MsgStream access and capabilities.
Definition: AthMessaging.h:55
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Preparation.mode
mode
Definition: Preparation.py:107
sTgcDigitMaker::GammaParameter::kParameter
double kParameter
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:77
sTgcDigitMaker::DigiConditions::rndEngine
CLHEP::HepRandomEngine * rndEngine
Definition: MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h:68
sTgcDigitMaker::addDigit
static void addDigit(sTgcDigitVec &digits, const Identifier &id, const uint16_t bctag, const double digittime, const double charge)
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:740
sTgcDigitMaker::DigiConditions::efficiencies
const Muon::DigitEffiData * efficiencies
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:60
compareGeometries.deltaX
float deltaX
Definition: compareGeometries.py:32
dot.dot
def dot(G, fn, nodesToHighlight=[])
Definition: dot.py:5
sTgcDigitMaker::Ionization::posOnWire
Amg::Vector3D posOnWire
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:87
sTgcDigitMaker::Ionization::posOnSegment
Amg::Vector3D posOnSegment
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:86
sTgcDigitMaker.h
tolerance
Definition: suep_shower.h:17
sTgcDigitMaker::m_doTimeOffsetStrip
bool m_doTimeOffsetStrip
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:155
PathResolver.h
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:28
sTgcDigitMaker::m_doPadSharing
bool m_doPadSharing
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:151
charge
double charge(const T &p)
Definition: AtlasPID.h:986
sTgcIdHelper
Definition: sTgcIdHelper.h:55
sTgcDigitMaker::getPadChargeFraction
static double getPadChargeFraction(double distance)
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:865
sTgcDigitMaker::Ionization::distance
double distance
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:85
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MuonGMR4::sTgcReadoutElement
Definition: MuonPhaseII/MuonDetDescr/MuonReadoutGeometryR4/MuonReadoutGeometryR4/sTgcReadoutElement.h:20
sTgcDigitMaker::DigiConditions::detMgr
const MuonGM::MuonDetectorManager * detMgr
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:59
MuonGMR4::sTgcReadoutElement::padDesign
const PadDesign & padDesign(const Identifier &measId) const
Retrieves the readoutElement Layer given the Identifier/Hash.
DeMoScan.index
string index
Definition: DeMoScan.py:362
python.SystemOfUnits.keV
float keV
Definition: SystemOfUnits.py:174
sTgcDigitMaker::m_idHelperSvc
const Muon::IMuonIdHelperSvc * m_idHelperSvc
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:139
sTgcDigitMaker::m_mostProbableArrivalTime
std::array< double, 5 > m_mostProbableArrivalTime
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:134
Muon::IMuonIdHelperSvc::toString
virtual std::string toString(const Identifier &id) const =0
print all fields to string
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path)
Definition: PathResolver.cxx:183
known
Definition: TrigBStoxAODTool.cxx:107
DeMoScan.first
bool first
Definition: DeMoScan.py:534
sTgcDigitMaker::m_theta
double m_theta
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:149
sTgcDigitMaker::digitMode
digitMode
Constructor initializing digitization parameters.
Definition: MuonPhaseII/MuonDigitization/sTgcDigitizationR4/sTgcDigitizationR4/sTgcDigitMaker.h:40
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
Muon::IMuonIdHelperSvc
Interface for Helper service that creates muon Identifiers and can be used to print Identifiers.
Definition: IMuonIdHelperSvc.h:27
sTgcDigitMaker::initialize
StatusCode initialize()
Initializes sTgcHitIdHelper and sTgcIdHelper, and call the functions to read files containing digitiz...
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:52
sTgcDigitMaker::m_StripResolution
double m_StripResolution
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:157
sTgcDigitMaker::sTgcDigitVec
std::vector< std::unique_ptr< sTgcDigit > > sTgcDigitVec
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:65
MuonGMR4::StripDesign::stripPitch
double stripPitch() const
Distance between two adjacent strips.
sTgcDigitMaker::sTgcDigitMaker
sTgcDigitMaker(const Muon::IMuonIdHelperSvc *idHelperSvc, const int channelTypes, double meanGasGain, bool doPadChargeSharing, double stripChargeScale, bool applyAsBuiltBLines)
Definition: MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx:34
sTgcDigitMaker::m_meanGasGain
double m_meanGasGain
Definition: MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h:150
MuonGMR4::PadDesign
Definition: PadDesign.h:24
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37
Identifier
Definition: IdentifierFieldParser.cxx:14
MuonGMR4::StripDesign::firstStripPos
const Amg::Vector2D & firstStripPos() const
Vector indicating the first strip position.