ATLAS Offline Software
Loading...
Searching...
No Matches
sTgcDigitMaker Class Reference

Handles the digitization of sTGC hits, converting simulated hits into sTGC digits. More...

#include <sTgcDigitMaker.h>

Inheritance diagram for sTgcDigitMaker:
Collaboration diagram for sTgcDigitMaker:

Classes

struct  DigiConditions
 Digitize a given hit, determining the time and charge spread on wires, pads and strips. More...
struct  GammaParameter
 Parameters of a gamma probability distribution function, required for estimating wire digit's time of arrival. More...
struct  Ionization
 Ionization object with distance, position on the hit segment and position on the wire. More...

Public Types

using sTgcDigitVec = std::vector<std::unique_ptr<sTgcDigit>>

Public Member Functions

 sTgcDigitMaker (const Muon::IMuonIdHelperSvc *idHelperSvc, const int channelTypes, double meanGasGain, bool doPadChargeSharing, bool applyAsBuiltBLines)
virtual ~sTgcDigitMaker ()
StatusCode initialize ()
 Initializes sTgcHitIdHelper and sTgcIdHelper, and call the functions to read files containing digitization parameters.
sTgcDigitVec executeDigi (const DigiConditions &condContainers, const sTGCSimHit &hit) const
bool msgLvl (const MSG::Level lvl) const
 Test the output level.
MsgStream & msg () const
 The standard message stream.
MsgStream & msg (const MSG::Level lvl) const
 The standard message stream.
void setLevel (MSG::Level lvl)
 Change the current logging level.

Private Member Functions

StatusCode readFileOfTimeArrival ()
 Read share/sTGC_Digitization_timeArrival.dat.
StatusCode readFileOfTimeOffsetStrip ()
 Read share/sTGC_Digitization_timeOffsetStrip.dat.
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.
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 cluster.
GammaParameter getGammaParameter (double distance) const
 Find the gamma pdf parameters of a given distance.
double getMostProbableArrivalTime (double distance) const
 Get the most probable time of arrival.
void initMessaging () const
 Initialize our message level and MessageSvc.

Static Private Member Functions

static void addDigit (sTgcDigitVec &digits, const Identifier &id, const uint16_t bctag, const double digittime, const double charge)
static double getPadChargeFraction (double distance)

Private Attributes

std::vector< GammaParameterm_gammaParameter
std::array< double, 5 > m_mostProbableArrivalTime {make_array<double, 5>(0.)}
std::array< double, 6 > m_timeOffsetStrip {make_array<double, 6>(0.)}
const Muon::IMuonIdHelperSvcm_idHelperSvc {nullptr}
int m_channelTypes {3}
 define offsets and widths of time windows for signals from wiregangs and strips.
double m_theta {10}
double m_meanGasGain {5.e4}
bool m_doPadSharing {false}
bool m_doTimeOffsetStrip {false}
double m_StripResolution {0.0949}
double m_posResIncident {1.}
double m_posResAngular {0.305/m_StripResolution}
double m_chargeAngularFactor {4.0}
bool m_applyAsBuiltBLines {false}
std::string m_nm
 Message source name.
boost::thread_specific_ptr< MsgStream > m_msg_tls
 MsgStream instance (a std::cout like with print-out levels)
std::atomic< IMessageSvc * > m_imsg { nullptr }
 MessageSvc pointer.
std::atomic< MSG::Level > m_lvl { MSG::NIL }
 Current logging level.
std::atomic_flag m_initialized ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
 Messaging initialized (initMessaging)

Static Private Attributes

static constexpr std::array< double, 2 > m_clusterProfile {0.573, 1.092}

Detailed Description

Handles the digitization of sTGC hits, converting simulated hits into sTGC digits.

for SCT_ChargeTrappingTool

All functionality of sTGC digitization is implemented to this class.

Definition at line 36 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

Member Typedef Documentation

◆ sTgcDigitVec

Constructor & Destructor Documentation

◆ sTgcDigitMaker()

sTgcDigitMaker::sTgcDigitMaker ( const Muon::IMuonIdHelperSvc * idHelperSvc,
const int channelTypes,
double meanGasGain,
bool doPadChargeSharing,
bool applyAsBuiltBLines )

Definition at line 33 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

38 : AthMessaging ("sTgcDigitMaker"),
39 m_idHelperSvc{idHelperSvc},
40 m_channelTypes{channelTypes},
41 m_meanGasGain{meanGasGain},
42 m_doPadSharing{doPadChargeSharing},
43 m_applyAsBuiltBLines(applyAsBuiltBLines) {}
AthMessaging()
Default constructor:
int m_channelTypes
define offsets and widths of time windows for signals from wiregangs and strips.

◆ ~sTgcDigitMaker()

sTgcDigitMaker::~sTgcDigitMaker ( )
virtualdefault

Member Function Documentation

◆ addDigit()

void sTgcDigitMaker::addDigit ( sTgcDigitVec & digits,
const Identifier & id,
const uint16_t bctag,
const double digittime,
const double charge )
staticprivate

Definition at line 735 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

739 {
740
741 constexpr double tolerance = 0.1;
742 if (std::find_if(digits.begin(),digits.end(), [&](std::unique_ptr<sTgcDigit>& known) {
743 return known->identify() == id && std::abs(digittime - known->time()) < tolerance;
744 }) == digits.end()) {
745 digits.push_back(std::make_unique<sTgcDigit>(id, bctag, digittime, charge, 0, 0));
746 }
747}
double charge(const T &p)
Definition AtlasPID.h:997
constexpr double tolerance

◆ executeDigi()

sTgcDigitMaker::sTgcDigitVec sTgcDigitMaker::executeDigi ( const DigiConditions & condContainers,
const sTGCSimHit & hit ) const

Definition at line 69 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

70 {
71 // SimHits without energy loss are not recorded.
72 double energyDeposit = hit.depositEnergy(); // Energy deposit in MeV
73 if(energyDeposit==0.) return {};
74
76 const sTgcIdHelper& idHelper{m_idHelperSvc->stgcIdHelper()};
77 sTgcSimIdToOfflineId simToOffline{&idHelper};
78 const Identifier offlineId = simToOffline.convert(hit.sTGCId());
79 const Identifier layid = idHelper.channelID(offlineId,
80 idHelper.multilayer(offlineId),
81 idHelper.gasGap(offlineId),
83
84 ATH_MSG_VERBOSE("sTgc hit: time " << hit.globalTime()
85 << " position " << Amg::toString(hit.globalPosition(), 2)
86 << " mclink " << hit.particleLink() << " PDG ID " << hit.particleEncoding() );
87
88
89 ATH_MSG_DEBUG("Retrieving detector element for: "<< m_idHelperSvc->toStringDetEl(layid) << " energyDeposit "<<energyDeposit );
90
91 const MuonGM::sTgcReadoutElement* detEl = condContainers.detMgr->getsTgcReadoutElement(layid);
92 if( !detEl ){ return {}; }
93
94
95 // DO THE DIGITIZATTION HERE ////////
96
97 // Required precision on length in mm
98 constexpr double length_tolerance = 0.01;
99
100 // Retrieve the wire surface
101 int surfHash_wire = detEl->surfaceHash(layid);
102 const Trk::PlaneSurface& SURF_WIRE = detEl->surface(surfHash_wire); // get the wire surface
103
104 // Hit post-Step global position
105 const Amg::Vector3D& GPOS{hit.globalPosition()};
106 // Hit pre-Step global position
107 const Amg::Vector3D& pre_pos{hit.globalPrePosition()};
108 // Hit global direction
109 const Amg::Vector3D& GLODIRE{hit.globalDirection()};
110
111 // Hit position in the wire surface's coordinate frame
112 Amg::Vector3D hitOnSurface_wire{SURF_WIRE.transform().inverse() * GPOS};
113 Amg::Vector3D pre_pos_wire_surf{SURF_WIRE.transform().inverse() * pre_pos};
114 Amg::Vector2D posOnSurf_wire{0.5 * (hitOnSurface_wire + pre_pos_wire_surf).block<2,1>(0,0)};
115
116 /* Determine the closest wire and the distance of closest approach
117 * Since most particles pass through the the wire plane between two wires,
118 * the nearest wire should be one of these two wire. Otherwise, the particle's
119 * trajectory is uncommon, and such rare case is not supported yet.
120 *
121 * Finding that nearest wire follows the following steps:
122 * - Compute the distance to the wire at the center of the current wire pitch
123 * - Compute the distance to the other adjacent wire
124 */
125
126 // Wire number of the current wire pitch
127 int wire_number = detEl->getDesign(layid)->wireNumber(posOnSurf_wire);
128
129 // If wire number is invalid, verify if hit is near the edge of the chamber.
130 const int number_wires = detEl->numberOfWires(layid);
131 if ((wire_number < 1) || (wire_number > number_wires)) {
132 // Compute the wire number using either the pos-step position or
133 // the pre-step position, whichever yields a valid wire number.
134 wire_number = detEl->getDesign(layid)->wireNumber(hitOnSurface_wire.block<2,1>(0,0));
135 if ((wire_number < 1) || (wire_number > number_wires)) {
136 wire_number = detEl->getDesign(layid)->wireNumber(pre_pos_wire_surf.block<2,1>(0,0));
137 }
138 }
139
140 // Compute the position of the ionization and its distance relative to a given sTGC wire.
141 Ionization ionization = pointClosestApproach(detEl, layid, wire_number, pre_pos_wire_surf, hitOnSurface_wire);
142 double dist_wire = ionization.distance;
143 if (dist_wire > 0.) {
144 // Determine the other adjacent wire, which is
145 // -1 if particle crosses the wire surface between wire_number-1 and wire_number
146 // +1 if particle crosses the wire surface between wire_number and wire_number+1
147 int adjacent = 1;
148 if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {adjacent = -1;}
149
150 // Find the position of the ionization with respect to the adjacent wire
151 Ionization ion_adj = pointClosestApproach(detEl, layid, wire_number+adjacent, pre_pos_wire_surf, hitOnSurface_wire);
152 // Keep the closest point
153 double dist_wire_adj = ion_adj.distance;
154 if ((dist_wire_adj > 0.) && (dist_wire_adj < dist_wire)) {
155 dist_wire = dist_wire_adj;
156 wire_number += adjacent;
157 ionization = std::move(ion_adj);
158 }
159 } else {
160 ATH_MSG_DEBUG("Failed to get the distance between the wire number = " << wire_number
161 << " and hit at " <<Amg::toString(hitOnSurface_wire, 2)
162 << ". Number of wires = " << number_wires
163 << ", "<<m_idHelperSvc->toStringGasGap(layid));
164 return {};
165 }
166
167 // Update the position of ionization on the wire surface
168 posOnSurf_wire = ionization.posOnWire.block<2,1>(0,0);
169 // Position of the ionization in the global coordinate frame
170 const Amg::Vector3D glob_ionization_pos = SURF_WIRE.transform() * ionization.posOnWire;
171
172 ATH_MSG_VERBOSE("Ionization_info: distance: " << ionization.distance
173 << " posOnTrack: " <<Amg::toString(ionization.posOnSegment, 3)
174 << " posOnWire: " << Amg::toString(ionization.posOnWire, 3)
175 << " hitGPos: " << Amg::toString(GPOS,3)
176 << " hitPrePos: " <<Amg::toString(pre_pos, 3)
177 << " EDep: " << hit.depositEnergy() << " EKin: " << hit.kineticEnergy()
178 << " pdgId: " << hit.particleEncoding()
179 <<" "<<m_idHelperSvc->toStringGasGap(layid));
180
181 // Distance should be in the range [0, 0.9] mm, excepting
182 // - particles pass through the wire plane near the edges
183 // - secondary particles created inside the gas gap that go through the gas gap partially.
184 // Most of such particles are not muons and have low kinetic energy.
185 // - particle with trajectory parallel to the sTGC wire plane
186 const double wire_pitch = detEl->wirePitch();
187 if ((dist_wire > 0.) && (std::abs(hit.particleEncoding()) == 13) && (dist_wire > (wire_pitch/2))) {
188 ATH_MSG_DEBUG("Distance to the nearest wire (" << dist_wire << ") is greater than expected.");
189 ATH_MSG_DEBUG("Hit globalPos: "<<Amg::toString(GPOS,2)
190 << " globalPrePos: " << Amg::toString(pre_pos, 2)
191 << " EDeposited: " << hit.depositEnergy()
192 << " EKinetic: " << hit.kineticEnergy()
193 << " pdgID: " << hit.particleEncoding()
194 << " "<<m_idHelperSvc->toStringGasGap(layid));
195 }
196
197 // Do not digitize hits that are too far from the nearest wire
198 if (dist_wire > wire_pitch) {
199 return {};
200 }
201
202 // Get the gamma pdf parameters associated with the distance of closest approach.
203 const GammaParameter gamParam = getGammaParameter(dist_wire);
204 const double par_kappa = gamParam.kParameter;
205 const double par_theta = gamParam.thetaParameter;
206 const double most_prob_time = getMostProbableArrivalTime(dist_wire);
207 // Compute the most probable value of the gamma pdf
208 double gamma_mpv = (par_kappa - 1) * par_theta;
209 // If the most probable value is less than zero, then set it to zero
210 if (gamma_mpv < 0.) {gamma_mpv = 0.;}
211 const double t0_par = most_prob_time - gamma_mpv;
212
213 // Digit time follows a gamma distribution, so a value val is
214 // chosen using a gamma random generator then is shifted by t0
215 // to account for drift time.
216 // Note: CLHEP::RandGamma takes the parameters k and lambda,
217 // where lambda = 1 / theta.
218 double digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
219
220 // Sometimes, digit_time is negative because t0_par can be negative.
221 // In such case, discard the negative value and shoot RandGamma for another value.
222 // However, if that has already been done many times then set digit_time to zero
223 // in order to avoid runaway loop.
224 constexpr int shoot_limit = 4;
225 int shoot_counter = 0;
226 while (digit_time < 0.) {
227 if (shoot_counter > shoot_limit) {
228 digit_time = 0.;
229 break;
230 }
231 digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
232 ++shoot_counter;
233 }
234
235 ATH_MSG_DEBUG("sTgcDigitMaker distance = " << dist_wire
236 << ", time = " << digit_time
237 << ", k parameter = " << par_kappa
238 << ", theta parameter = " << par_theta
239 << ", most probable time = " << most_prob_time);
240
241 bool isValid = false;
243 if (condContainers.efficiencies) {
244 const double efficiency = condContainers.efficiencies->getEfficiency(layid);
245 // Lose Hits to match HV efficiency
246 if (CLHEP::RandFlat::shoot(condContainers.rndmEngine,0.0,1.0) > efficiency) return {};
247 }
248
249
250 sTgcDigitVec digits{};
251
252 //const double stripPropagationTime = 3.3*CLHEP::ns/CLHEP::m * detEl->distanceToReadout(posOnSurf_strip, elemId); // 8.5*ns/m was used until MC10.
253 const double stripPropagationTime = 0.; // 8.5*ns/m was used until MC10.
254
255 double sDigitTimeWire = digit_time;
256 double sDigitTimePad = sDigitTimeWire;
257 double sDigitTimeStrip = sDigitTimeWire + stripPropagationTime;
258
259 uint16_t bctag = 0;
260
261 //##################################################################################
262 //######################################### strip readout ##########################
263 //##################################################################################
264 ATH_MSG_DEBUG("sTgcDigitMaker::strip response ");
266
267 Identifier newId = idHelper.channelID(layid,idHelper.multilayer(layid),
268 idHelper.gasGap(layid), channelType, 1, isValid);
269
270 // get strip surface
271 int surfHash_strip = detEl->surfaceHash(newId);
272 const Trk::PlaneSurface& SURF_STRIP = detEl->surface(surfHash_strip); // get the strip surface
273
274 const Amg::Vector3D hitOnSurface_strip = SURF_STRIP.transform().inverse()*glob_ionization_pos;
275 Amg::Vector2D posOnSurf_strip {Amg::Vector2D::Zero()};
276
278 //This block is used to apply As-Built and BLine corrections for dedicated studies.
279 Amg::Vector3D posAfterAsBuilt {Amg::Vector3D::Zero()};
280 detEl->spacePointPosition(newId, hitOnSurface_strip.x(), hitOnSurface_strip.y(), posAfterAsBuilt);
281 posOnSurf_strip = posAfterAsBuilt.block<2,1>(0,0);
282 } else {
283 posOnSurf_strip = hitOnSurface_strip.block<2,1>(0,0);
284 }
285
286 bool insideBounds = SURF_STRIP.insideBounds(posOnSurf_strip);
287 if(!insideBounds) {
288 ATH_MSG_DEBUG("Outside of the strip surface boundary : " << m_idHelperSvc->toString(newId) << "; local position " <<posOnSurf_strip );
289 return {};
290 }
291
292 //************************************ find the nearest readout element **************************************
293
294 int stripNumber = detEl->stripNumber(posOnSurf_strip, newId);
295 if( stripNumber == -1 ){
296 // Verify if the energy deposit is at the boundary
297 const double new_posX = (posOnSurf_strip.x() > 0.0)? posOnSurf_strip.x() - length_tolerance
298 : posOnSurf_strip.x() + length_tolerance;
299 const Amg::Vector2D new_position(new_posX, posOnSurf_strip.y());
300 stripNumber = detEl->stripNumber(new_position, newId);
301 // Skip hit if still unable to obtain strip number
302 if (stripNumber < 1) {
303 ATH_MSG_WARNING("Failed to obtain strip number " << m_idHelperSvc->toString(newId) );
304 ATH_MSG_WARNING("Position on strip surface = (" << posOnSurf_strip.x() << ", " << posOnSurf_strip.y() << ")");
305 return {};
306 }
307 }
308 newId = idHelper.channelID(layid, idHelper.multilayer(layid),
309 idHelper.gasGap(layid), channelType, stripNumber, isValid);
310 if(!isValid && stripNumber != -1) {
311 ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(newId) );
312 return {};
313 }
314
315 const int NumberOfStrips = detEl->numberOfStrips(newId);
316 const double stripHalfPitch = detEl->channelPitch(newId)*0.5; // 3.2/2 = 1.6 mm
317
318 //************************************ conversion of energy to charge **************************************
319
320 // Typical ionized charge in pC per keV deposited. The constant is determined from ionization
321 // study with Garfield program. A note titled "Charge Energy Relation" which outlines
322 // conversion can be found here:
323 // https://cernbox.cern.ch/index.php/apps/files/?dir=/__myshares/Documents (id:274113) // link is dead
324 const double ionized_charge = (5.65E-6)*energyDeposit/CLHEP::keV;
325
326 // To get avalanche gain, polya function is taken from Blum paper https://inspirehep.net/literature/807304
327 // 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);
328
329 // Mean value for total gain due to E field;
330 // To calculate this gain from polya distibution, we replace in gamma PDF:
331 // alpha = 1+theta and
332 // beta = 1+theta/mean
333 // With these substitutions, gamma PDF gives the same sampling values as those from polya PDF.
334 const double gain = CLHEP::RandGamma::shoot(condContainers.rndmEngine, 1. + m_theta, (1. + m_theta)/m_meanGasGain);
335
336 // total charge after avalanche
337 const double total_charge = gain*ionized_charge;
338
339 //************************************ spread charge among readout element **************************************
340 // Charge Spread including tan(theta) resolution term.
341 const double tan_theta = GLODIRE.perp()/GLODIRE.z();
342 // The angle dependance on strip resolution goes as tan^2(angle)
343 const double angle_dependency = std::hypot(m_posResIncident, m_posResAngular * tan_theta);
344
345 const double cluster_posX = posOnSurf_strip.x();
346 double peak_position = CLHEP::RandGaussZiggurat::shoot(condContainers.rndmEngine, cluster_posX, m_StripResolution*angle_dependency);
347
348 // Each readout plane reads about half the total charge produced on the wire,
349 // including a tan(theta) term to describe the increase of charge with incident angle
350 const double norm = 0.5 * total_charge;
351 // Strip cluster charge profile described by a double Gaussian.
352 // The independent parameter x is in the units of strip channel, one strip channel = 3.2 mm,
353 // so convert position from mm to strip channel if it is not already.
354std::unique_ptr<TF1> clusterProfile = std::make_unique<TF1>("fgaus",
355 "0.5 * [2] / (sqrt(2 * TMath::Pi()) * [0]) * exp(-0.5 * (x / [0])^2) + "
356 "0.5 * [2] / (sqrt(2 * TMath::Pi()) * [1]) * exp(-0.5 * (x / [1])^2)", -300., 300.);
357 clusterProfile->SetParameters(m_clusterProfile[0], // sigma of 1st Gaussian
358 m_clusterProfile[1], // sigma of 2nd Gaussian
359 norm); // normalization factor
360
361 // Lower limit on strip charge (arbitrary limit), in pC, which has the same units as the parameter ionized_charge.
362 constexpr double tolerance_charge = 0.0005;
363 // Set a maximum number of neighbour strips to avoid very long loop. This is an arbitrary number.
364 constexpr unsigned int max_neighbor = 10;
365
366 // Spread charge on the strips that are on the upper half of the strip cluster
367 for (unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
368 int currentStrip = stripNumber + iStrip;
369 if (currentStrip > NumberOfStrips) break;
370
371 // Get the strip identifier and create the digit
372 newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid),
373 channelType, currentStrip, isValid);
374 if (isValid) {
375 Amg::Vector2D locpos{Amg::Vector2D::Zero()};
376 if (!detEl->stripPosition(newId, locpos)) {
377 ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(newId) );
378 }
379
380 // Estimate the digit charge
381 // Position with respect to the peak of the charge curve
382 double x_relative = locpos.x() - peak_position;
383 // In clusterProfile curve, position should be in the units of strip channel
384 double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
385 // If charge is too small, stop creating neighbor strip
386 if (charge < tolerance_charge) break;
387
388 // Estimate digit time
389 double strip_time = sDigitTimeStrip;
390 // Strip time response can be delayed due to the resistive layer.
391 // A correction would be required if the actual VMM front-end doesn't re-align the strip timing.
393 // Determine how far the current strip is from the middle strip
394 int indexFromMiddleStrip = iStrip;
395 if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
396 indexFromMiddleStrip = iStrip - 1;
397 // Add time delay due to resistive layer
398 strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
399 }
400
401 addDigit(digits, newId, bctag, strip_time, charge);
402
403 ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
404 << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
405 << ", neighbor index = " << iStrip
406 << ", strip position = "<<Amg::toString(locpos, 2));
407 }
408 }
409
410 // The lower half of the strip cluster
411 for (unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
412 int currentStrip = stripNumber - iStrip;
413 if (currentStrip < 1) break;
414
415 newId = idHelper.channelID(layid,idHelper.multilayer(layid),
416 idHelper.gasGap(layid), channelType, currentStrip, isValid);
417 if (isValid) {
418 Amg::Vector2D locpos{Amg::Vector2D::Zero()};
419 if (!detEl->stripPosition(newId, locpos)) {
420 ATH_MSG_WARNING("Failed to obtain local position for identifier " << m_idHelperSvc->toString(newId) );
421 }
422
423 // Estimate the digit charge
424 double x_relative = locpos.x() - peak_position;
425 double charge = std::hypot(1, m_chargeAngularFactor * tan_theta) * clusterProfile->Integral(x_relative/(2*stripHalfPitch) - 0.5, x_relative/(2*stripHalfPitch) + 0.5);
426 if (charge < tolerance_charge) break;
427
428 // Estimate digit time
429 double strip_time = sDigitTimeStrip;
430 // Time delay due to resistive layer
432 int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
433 strip_time += getTimeOffsetStrip(indexFromMiddleStrip);
434 }
435
436 addDigit(digits, newId, bctag, strip_time, charge);
437
438 ATH_MSG_VERBOSE("Created a strip digit: strip number = " << currentStrip << ", charge = " << charge
439 << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
440 << ", neighbor index = " << iStrip
441 << ", strip position = " << Amg::toString(locpos, 2));
442 }
443 }
444 // end of strip digitization
445
446 if(m_channelTypes==1) {
447 ATH_MSG_WARNING("Only digitize strip response !");
448 return digits;
449 }
450
451 //##################################################################################
452 //######################################### pad readout ##########################
453 //##################################################################################
454 ATH_MSG_DEBUG("sTgcDigitMaker::pad response ");
456
457 Identifier PAD_ID = idHelper.channelID(layid, idHelper.multilayer(layid),
458 idHelper.gasGap(layid), channelType, 1);// find the a pad id
459 //************************************ find the nearest readout element **************************************
460 int surfHash_pad = detEl->surfaceHash(PAD_ID);
461 const Trk::PlaneSurface& SURF_PAD = detEl->surface(surfHash_pad); // get the pad surface
462
463 const Amg::Vector3D hitOnSurface_pad = SURF_PAD.transform().inverse()*glob_ionization_pos;
464 const Amg::Vector2D posOnSurf_pad(hitOnSurface_pad.x(), hitOnSurface_pad.y());
465
466
467 insideBounds = SURF_PAD.insideBounds(posOnSurf_pad);
468
469 if(insideBounds) {
470 int padNumber = detEl->stripNumber(posOnSurf_pad, PAD_ID);
471 if( padNumber == -1 ){
472 // Verify if the energy deposit is at the boundary
473 const double new_posX = (posOnSurf_pad.x()>0.0)? posOnSurf_pad.x()-length_tolerance
474 : posOnSurf_pad.x()+length_tolerance;
475 const double new_posY = (posOnSurf_pad.y()>0.0)? posOnSurf_pad.y()-length_tolerance
476 : posOnSurf_pad.y()+length_tolerance;
477 const Amg::Vector2D new_position(new_posX, new_posY);
478 padNumber = detEl->stripNumber(new_position, PAD_ID);
479 // Skip hit if still unable to obtain pad number
480 if (padNumber < 1) {
481 ATH_MSG_WARNING("Failed to obtain pad number " << m_idHelperSvc->toString(PAD_ID) );
482 ATH_MSG_WARNING("Position on pad surface = (" << posOnSurf_pad.x() << ", " << posOnSurf_pad.y() << ")");
483 return digits;
484 }
485 }
486 newId = idHelper.channelID(layid, idHelper.multilayer(layid),
487 idHelper.gasGap(layid), channelType, padNumber, isValid);
488 if(isValid) {
489 // Find centre position of pad
490 Amg::Vector2D padPos{Amg::Vector2D::Zero()};
491 if (!detEl->stripPosition(newId, padPos)) {
492 ATH_MSG_ERROR("Failed to obtain local position for neighbor pad with identifier " << m_idHelperSvc->toString(newId) );
493 return digits;
494 }
495 // Pad sharing needs to look at position on hit vs pad boundaries
496 const Amg::Vector2D diff = posOnSurf_pad - padPos;
497 double halfPadWidthX = 0.5*detEl->getPadDesign(newId)->channelWidth(posOnSurf_pad, true, true);
498 double halfPadWidthY = 0.5*detEl->getPadDesign(newId)->channelWidth(posOnSurf_pad, false);
499
500 // Charge sharing happens within 4mm window for pads
501 // i.e. the charge is spread over a 2D gaussian of total radius about 4mm
502 // This value depends on actual width of the distribution, but we don't have the
503 // width of the charge distribution for pads. So lets consider 2.5*sigma, where
504 // sigma is the width of the charge distribution for strips.
505 double deltaX = halfPadWidthX - std::abs(diff.x());
506 double deltaY= halfPadWidthY - std::abs(diff.y());
507 bool isNeighX = deltaX < 2.5*m_clusterProfile[1];
508 bool isNeighY = deltaY < 2.5*m_clusterProfile[1];
509 // Pad width can be calculated to be very slightly larger than it should due to rounding errors
510 // 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.
511 if (deltaX < 0.) deltaX = 0.1;
512 if (deltaY < 0.) deltaY = 0.1;
513
514 if (m_doPadSharing && (isNeighX || isNeighY)){
515 std::pair<int, int> padEtaPhi(detEl->getPadDesign(newId)->channelNumber(padPos));
516 // Phi == 1 at the right in the local geometry
517 // In local coordinates, the pad to the right has phi' = phi - 1
518 int newPhi = diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
519 int newEta = diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
520 bool validEta = newEta > 0 && newEta < detEl->getPadDesign(newId)->nPadH+1;
521 bool validPhi = newPhi > 0 && newPhi < detEl->getPadDesign(newId)->nPadColumns+1;
522
523 if (isNeighX && isNeighY && validEta && validPhi){
524 // 4 pads total; makes life a bit harder. Corner of 4 valid pads
525 Identifier neigh_ID_X = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
526 channelType, padEtaPhi.first, newPhi, isValid);
527 Identifier neigh_ID_Y = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
528 channelType, newEta, padEtaPhi.second, isValid);
529 Identifier neigh_ID_XY = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
530 channelType, newEta, newPhi, isValid);
531 double xQfraction = getPadChargeFraction(deltaX);
532 double yQfraction = getPadChargeFraction(deltaY);
533
534 // Main pad gets 1 - Qfraction of total
535 addDigit(digits, neigh_ID_X, bctag, sDigitTimePad, xQfraction*(1.-yQfraction)*0.5*total_charge);
536 addDigit(digits, neigh_ID_Y, bctag, sDigitTimePad, yQfraction*(1.-xQfraction)*0.5*total_charge);
537 addDigit(digits, neigh_ID_XY, bctag, sDigitTimePad, xQfraction*yQfraction*0.5*total_charge);
538 addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction-yQfraction+xQfraction*yQfraction)*0.5*total_charge);
539 } else if (isNeighX && validPhi){
540 // There is only 1 neighbor, immediately to the left or right.
541 Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
542 channelType, padEtaPhi.first, newPhi, isValid);
543
544 // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
545 // Main pad gets 1 - Qfraction of total
546 double xQfraction = getPadChargeFraction(deltaX);
547 addDigit(digits, newId, bctag, sDigitTimePad, (1.-xQfraction)*0.5*total_charge);
548 addDigit(digits, neigh_ID, bctag, sDigitTimePad, xQfraction*0.5*total_charge);
549 }
550 else if (isNeighY && validEta){
551 // There is only 1 neighbor, immediately above or below
552 Identifier neigh_ID = idHelper.padID(newId, idHelper.multilayer(newId), idHelper.gasGap(newId),
553 channelType, newEta, padEtaPhi.second, isValid);
554
555 // NeighborPad gets Qfraction of the total pad charge: 0.5*total_charge
556 // Main pad gets 1 - Qfraction of total
557 double yQfraction = getPadChargeFraction(deltaX);
558 addDigit(digits, newId, bctag, sDigitTimePad, (1.-yQfraction)*0.5*total_charge);
559 addDigit(digits, neigh_ID, bctag, sDigitTimePad, yQfraction*0.5*total_charge);
560 }
561
562 }
563 else{
564 // No charge sharing: hit is nicely isolated within pad
565 addDigit(digits, newId, bctag, sDigitTimePad, 0.5*total_charge);
566 }
567
568 }
569 else if(padNumber != -1) {
570 ATH_MSG_ERROR("Failed to obtain identifier " << m_idHelperSvc->toString(newId) );
571 }
572 }
573 else {
574 ATH_MSG_DEBUG("Outside of the pad surface boundary :" << m_idHelperSvc->toString(PAD_ID)<< " local position " <<posOnSurf_pad );
575 }
576
577 if(m_channelTypes==2) {
578 ATH_MSG_WARNING("Only digitize strip/pad response !");
579 return digits;
580 }
581
582
583 //##################################################################################
584 //######################################### wire readout ##########################
585 //##################################################################################
586 ATH_MSG_DEBUG("sTgcDigitMaker::wire response ");
588
589 // Find the ID of the first wiregroup
590 Identifier WIREGP_ID = idHelper.channelID(layid, idHelper.multilayer(layid),
591 idHelper.gasGap(layid), channelType, 1);
592
593 //************************************ find the nearest readout element **************************************
594 insideBounds = SURF_WIRE.insideBounds(posOnSurf_wire);
595
596 if(insideBounds) {
597 // Determine the wire number
598
599 int wiregroupNumber = detEl->stripNumber(posOnSurf_wire, WIREGP_ID);
600 if( wiregroupNumber == -1 ){
601 // Verify if the energy deposit is at the boundary
602 const double new_posX = (posOnSurf_wire.x() > 0.0)? posOnSurf_wire.x() - length_tolerance
603 : posOnSurf_wire.x() + length_tolerance;
604 const Amg::Vector2D new_position(new_posX, posOnSurf_wire.y());
605 wiregroupNumber = detEl->stripNumber(new_position, WIREGP_ID);
606 // Skip hit if still unable to obtain pad number
607 if (wiregroupNumber < 1) {
608 ATH_MSG_WARNING("Failed to obtain wire number " << m_idHelperSvc->toString(WIREGP_ID) );
609 ATH_MSG_WARNING("Position on wire surface = (" << posOnSurf_wire.x() << ", " << posOnSurf_wire.y() << ")");
610 return digits;
611 }
612 }
613
614 // Find ID of the actual wiregroup
615 newId = idHelper.channelID(layid, idHelper.multilayer(layid), idHelper.gasGap(layid),
616 channelType, wiregroupNumber, isValid);
617
618 if(isValid) {
619 int NumberOfWiregroups = detEl->numberOfStrips(newId);
620 if(wiregroupNumber>=1&&wiregroupNumber<=NumberOfWiregroups) {
621 addDigit(digits, newId, bctag, sDigitTimeWire, total_charge);
622 }
623 } // end of if(isValid)
624 else if (wiregroupNumber != -1){
625 ATH_MSG_ERROR("Failed to obtain wiregroup identifier " << m_idHelperSvc->toString(newId) );
626 }
627 }
628 else {
629 ATH_MSG_DEBUG("Outside of the wire surface boundary :" << m_idHelperSvc->toString(WIREGP_ID)<< " local position " <<posOnSurf_wire );
630 }
631 // end of wire digitization
632
633 return digits;
634}
#define ATH_MSG_ERROR(x)
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
Definition AtlasPID.h:878
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
virtual const Trk::PlaneSurface & surface() const override
access to chamber surface (phi orientation), uses the first gas gap
const MuonPadDesign * getPadDesign(const Identifier &id) const
returns the MuonChannelDesign class for the given identifier
virtual bool stripPosition(const Identifier &id, Amg::Vector2D &pos) const override final
strip position - should be renamed to channel position If the strip number is outside the range of va...
const MuonChannelDesign * getDesign(const Identifier &id) const
returns the MuonChannelDesign class for the given identifier
virtual int numberOfStrips(const Identifier &layerId) const override final
number of strips per layer
int numberOfWires(const Identifier &id) const
Get the total number of wires (single wires) of a chamber.
double channelPitch(const Identifier &id) const
Channel pitch.
virtual int surfaceHash(const Identifier &id) const override final
returns the hash to be used to look up the surface and transform in the MuonClusterReadoutElement tra...
virtual bool spacePointPosition(const Identifier &phiId, const Identifier &etaId, Amg::Vector2D &pos) const override final
space point position for a given pair of phi and eta identifiers The LocalPosition is expressed in th...
virtual int stripNumber(const Amg::Vector2D &pos, const Identifier &id) const override final
strip number corresponding to local position.
virtual bool insideBounds(const Amg::Vector2D &locpos, double tol1=0., double tol2=0.) const override
This method calls the inside() method of the Bounds.
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
const Amg::Vector3D & globalPrePosition() const
Definition sTGCSimHit.h:49
HitID sTGCId() const
Definition sTGCSimHit.h:51
double globalTime() const
Definition sTGCSimHit.h:41
const Amg::Vector3D & globalDirection() const
Definition sTGCSimHit.h:46
double kineticEnergy() const
Definition sTGCSimHit.h:48
const Amg::Vector3D & globalPosition() const
Definition sTGCSimHit.h:44
const HepMcParticleLink & particleLink() const
Definition sTGCSimHit.h:96
int particleEncoding() const
Definition sTGCSimHit.h:45
double depositEnergy() const
Definition sTGCSimHit.h:47
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...
std::vector< std::unique_ptr< sTgcDigit > > sTgcDigitVec
GammaParameter getGammaParameter(double distance) const
Find the gamma pdf parameters of a given distance.
static void addDigit(sTgcDigitVec &digits, const Identifier &id, const uint16_t bctag, const double digittime, const double charge)
double getMostProbableArrivalTime(double distance) const
Get the most probable time of arrival.
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.
int multilayer(const Identifier &id) const
Identifier padID(int stationName, int stationEta, int stationPhi, int multilayer, int gasGap, int channelType, int padEta, int padPhi) const
int gasGap(const Identifier &id) const override
get the hashes
Identifier channelID(int stationName, int stationEta, int stationPhi, int multilayer, int gasGap, int channelType, int channel) const
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="")
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
bool adjacent(unsigned int strip1, unsigned int strip2)
@ energyDeposit
setWord1 uint16_t
int wireNumber(const Amg::Vector2D &pos) const
calculate the sTGC wire number. The method can return a value outside the range [1,...
std::pair< int, int > channelNumber(const Amg::Vector2D &pos) const
calculate local channel number, range 1=nstrips like identifiers.
double channelWidth(const Amg::Vector2D &pos, bool measPhi, bool preciseMeas=false) const
calculate local channel width
Parameters of a gamma probability distribution function, required for estimating wire digit's time of...
Ionization object with distance, position on the hit segment and position on the wire.
Identifier convert(int simId) const

◆ getGammaParameter()

sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter ( double distance) const
private

Find the gamma pdf parameters of a given distance.

Definition at line 790 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

790 {
791
792 const double d = std::abs(distance);
793 // Find the parameters assuming the container is sorted in ascending order of 'lowEdge'
794 int index{-1};
795 for (const auto& par: m_gammaParameter) {
796 if (d < par.lowEdge) {
797 break;
798 }
799 ++index;
800 }
801 return m_gammaParameter.at(index);
802}
str index
Definition DeMoScan.py:362

◆ getMostProbableArrivalTime()

double sTgcDigitMaker::getMostProbableArrivalTime ( double distance) const
private

Get the most probable time of arrival.

Definition at line 805 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

805 {
806
807 const double d{std::abs(distance)};
808 double mpt{0};
809 for (size_t t = 0 ; t < m_mostProbableArrivalTime.size(); ++t){
810 mpt += m_mostProbableArrivalTime[t] * std::pow(d, t);
811 }
812 return mpt;
813}

◆ getPadChargeFraction()

double sTgcDigitMaker::getPadChargeFraction ( double distance)
staticprivate

Definition at line 860 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

860 {
861 // The charge fraction that is found past a distance x away from the
862 // centre of a 2D gaussian distribution of width of cluster profile is
863 // described by a modified error function.
864
865 // The modified error function perfectly describes
866 // the pad charge sharing distribution figure 16 of the sTGC
867 // testbeam paper https://arxiv.org/pdf/1509.06329.pdf
868
869 return 0.5 * (1.0 - std::erf( distance / (std::sqrt(2) * m_clusterProfile[1])));
870}

◆ getTimeOffsetStrip()

double sTgcDigitMaker::getTimeOffsetStrip ( size_t neighbor_index) const
private

Get digit time offset of a strip depending on its relative position to the strip at the centre of the cluster.

It returns 0 ns by default, as well as when it fails or container is empty.

Definition at line 855 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

855 {
856 return m_timeOffsetStrip.at(std::min(neighbor_index, m_timeOffsetStrip.size() -1));
857}

◆ initialize()

StatusCode sTgcDigitMaker::initialize ( )

Initializes sTgcHitIdHelper and sTgcIdHelper, and call the functions to read files containing digitization parameters.

Definition at line 49 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

49 {
50
51 // Read share/sTGC_Digitization_timeArrivale.dat, containing the digit time of arrival
53
54 // Read share/sTGC_Digitization_timeOffsetStrip.dat if the the strip time correction is enable
57 }
58 // check the digitization channel type
60 ATH_MSG_ERROR("Invalid ChannelTypes : " << m_channelTypes << " (valid values are : 1 --> strips ; 2 --> strips / wires ; 3 --> strips / wires / pads)");
61 return StatusCode::FAILURE;
62 }
63 return StatusCode::SUCCESS;
64}
#define ATH_CHECK
Evaluate an expression and check for errors.
StatusCode readFileOfTimeOffsetStrip()
Read share/sTGC_Digitization_timeOffsetStrip.dat.
StatusCode readFileOfTimeArrival()
Read share/sTGC_Digitization_timeArrival.dat.

◆ initMessaging()

void AthMessaging::initMessaging ( ) const
privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

40{
42 // If user did not set an explicit level, set a default
43 if (m_lvl == MSG::NIL) {
44 m_lvl = m_imsg ?
45 static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
46 MSG::INFO;
47 }
48}
std::string m_nm
Message source name.
std::atomic< IMessageSvc * > m_imsg
MessageSvc pointer.
std::atomic< MSG::Level > m_lvl
Current logging level.
IMessageSvc * getMessageSvc(bool quiet=false)

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 163 of file AthMessaging.h.

164{
165 MsgStream* ms = m_msg_tls.get();
166 if (!ms) {
167 if (!m_initialized.test_and_set()) initMessaging();
168 ms = new MsgStream(m_imsg,m_nm);
169 m_msg_tls.reset( ms );
170 }
171
172 ms->setLevel (m_lvl);
173 return *ms;
174}
boost::thread_specific_ptr< MsgStream > m_msg_tls
MsgStream instance (a std::cout like with print-out levels)
void initMessaging() const
Initialize our message level and MessageSvc.

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 178 of file AthMessaging.h.

179{ return msg() << lvl; }
MsgStream & msg() const
The standard message stream.

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl) const
inlineinherited

Test the output level.

Parameters
lvlThe message level to test against
Returns
boolean Indicating if messages at given level will be printed
Return values
trueMessages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

152{
153 if (m_lvl <= lvl) {
154 msg() << lvl;
155 return true;
156 } else {
157 return false;
158 }
159}

◆ pointClosestApproach()

sTgcDigitMaker::Ionization sTgcDigitMaker::pointClosestApproach ( const MuonGM::sTgcReadoutElement * readoutEle,
const Identifier & id,
int wireNumber,
const Amg::Vector3D & preStepPos,
const Amg::Vector3D & postStepPos ) const
private

Determine the points where the distance between two segments is smallest.

Given two segments, e.g. a particle trajectory and a sTGC wire, solve for the two points, the point on the trajectory and the point on the wire, where the distance between the two segments is the smallest.

Positions returned are in the local coordinate frame of the wire plane. Returns an object with distance of -9.99 in case of error.

Definition at line 637 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

641 {
642 // tolerance in mm
643 constexpr double tolerance = 0.001;
644 // minimum segment length
645 constexpr double min_length = 0.1;
646
647 // Position of the ionization
648 Ionization ionization;
649
650 // Wire number should be one or greater
651 if (wireNumber < 1) return ionization;
652
653 // Wire number too large
654 if (wireNumber > detEl->numberOfWires(id)) return ionization;
655
656 // Finding smallest distance and the points at the smallest distance.
657 // The smallest distance between two lines is perpendicular to both lines.
658 // We can construct two lines in the wire surface local coordinate frame:
659 // - one for the hit segment with equation h0 + t * v_h, where h0 is a point
660 // and v_h is the unit vector of the hit segment
661 // - another for the wire with similar equation w0 + s * v_w, where w0 is a
662 // point and v_w is the unit vector of the wire line
663 // Then it is possible to determine the factors s and t by requiring the
664 // dot product to be zero:
665 // 1. (h0 - w0) \dot v_h = 0
666 // 2. (h0 - w0) \dot v_w = 0
667
668 // Wire pitch
669 const double wire_pitch = detEl->wirePitch();
670 // Wire local position on the wire plane, the y-coordinate is arbitrary and z-coordinate is zero
671 const double wire_posX = detEl->positionFirstWire(id) + (wireNumber - 1) * wire_pitch;
672 const Amg::Vector3D wire_position(wire_posX, postStepPos.y(), 0.);
673 // The wires are parallel to Y in the wire plane's coordinate frame
674 const Amg::Vector3D wire_direction{Amg::Vector3D::UnitY()};
675
676 // particle trajectory
677 Amg::Vector3D hit_direction{postStepPos - preStepPos};
678 const double seg_length = hit_direction.mag();
679 if (seg_length > tolerance) hit_direction /= seg_length;
680
681 // Find the point on the track segment that is closest to the wire
682 if (seg_length < min_length) {
683 ionization.posOnSegment = postStepPos;
684 ionization.posOnWire = wire_position;
685 ionization.distance = std::hypot(postStepPos.x() - wire_posX, postStepPos.z());
686 return ionization;
687 }
688
689 // Dot product between the wire and hit direction
690 const double cos_theta = wire_direction.dot(hit_direction);
691 // distance between the hit and wire
692 const Amg::Vector3D dist_wire_hit = postStepPos - wire_position;
693
694 // Verifier the special case where the two lines are parallel
695 if (std::abs(cos_theta - 1.0) < tolerance) {
696 ATH_MSG_DEBUG("The track segment is parallel to the wire, position of digit is undefined");
697 ionization.posOnSegment = 0.5 * (postStepPos + preStepPos);
698 ionization.posOnWire = Amg::Vector3D(wire_posX, ionization.posOnSegment.y(), 0.0);
699 ionization.distance = std::hypot(ionization.posOnSegment.x() - wire_posX,
700 ionization.posOnSegment.z());
701 return ionization;
702 }
703
704 // Perpendicular component squared
705 const double sin_theta_2 = 1.0 - cos_theta * cos_theta;
706
707 //* Point on the hit segment
708 const double factor_hit = (-dist_wire_hit.dot(hit_direction)
709 + dist_wire_hit.dot(wire_direction) * cos_theta) / sin_theta_2;
710 Amg::Vector3D ionization_pos = postStepPos + factor_hit * hit_direction;
711
712 // If the point is on the track segment, then compute the other point on the wire.
713 // Otherwise, set the ionization at the pre-step position and compute where it
714 // should be on the wire.
715 Amg::Vector3D pos_on_wire{Amg::Vector3D::Zero()};
716 if (factor_hit < seg_length) {
717 //* Point on the wire line
718 const double factor_wire = (dist_wire_hit.dot(wire_direction)
719 - dist_wire_hit.dot(hit_direction) * cos_theta) / sin_theta_2;
720 pos_on_wire = wire_position + factor_wire * wire_direction;
721 } else {
722 ionization_pos = preStepPos;
723 pos_on_wire = wire_position - (seg_length * cos_theta) * wire_direction;
724 }
725
726 // Save the distance and ionization position
727 ionization.distance = (ionization_pos - pos_on_wire).mag();
728 ionization.posOnSegment = ionization_pos;
729 ionization.posOnWire = pos_on_wire;
730
731 return ionization;
732}
Scalar mag() const
mag method

◆ readFileOfTimeArrival()

StatusCode sTgcDigitMaker::readFileOfTimeArrival ( )
private

Read share/sTGC_Digitization_timeArrival.dat.

Definition at line 750 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

750 {
751 // Verify the file sTGC_Digitization_timeArrival.dat exists
752 const std::string file_name = "sTGC_Digitization_timeArrival.dat";
753 std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
754 if(file_path.empty()) {
755 ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
756 return StatusCode::FAILURE;
757 }
758
759 // Open the sTGC_Digitization_timeArrival.dat file
760 std::ifstream ifs{file_path, std::ios::in};
761 if(ifs.bad()) {
762 ATH_MSG_FATAL("sTgcDigitMaker: Failed to open time of arrival file " << file_name );
763 return StatusCode::FAILURE;
764 }
765
766 // Read the sTGC_Digitization_timeWindowOffset.dat file
767 std::string line;
768 GammaParameter param{};
769
770 while (std::getline(ifs, line)) {
771 std::string key;
772 std::istringstream iss(line);
773 iss >> key;
774 if (key.compare("bin") == 0) {
775 iss >> param.lowEdge >> param.kParameter >> param.thetaParameter;
776 m_gammaParameter.push_back(param);
777 } else if (key.compare("mpv") == 0) {
778 double mpt;
779 int idx{0};
780 while (iss >> mpt) {m_mostProbableArrivalTime[idx++] = mpt;}
781 }
782 }
783
784 // Close the file
785 ifs.close();
786 return StatusCode::SUCCESS;
787}
#define ATH_MSG_FATAL(x)
static std::string find_file(const std::string &logical_file_name, const std::string &search_path)
file_path
Definition athena.py:94

◆ readFileOfTimeOffsetStrip()

StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip ( )
private

Read share/sTGC_Digitization_timeOffsetStrip.dat.

Definition at line 816 of file MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx.

816 {
817 // Verify the file sTGC_Digitization_timeOffsetStrip.dat exists
818 const std::string file_name = "sTGC_Digitization_timeOffsetStrip.dat";
819 std::string file_path = PathResolver::find_file(file_name, "DATAPATH");
820 if(file_path.empty()) {
821 ATH_MSG_FATAL("readFileOfTimeWindowOffset(): Could not find file " << file_name );
822 return StatusCode::FAILURE;
823 }
824
825 // Open the sTGC_Digitization_timeOffsetStrip.dat file
826 std::ifstream ifs{file_path, std::ios::in};
827 if(ifs.bad()) {
828 ATH_MSG_FATAL("Failed to open time of arrival file " << file_name );
829 return StatusCode::FAILURE;
830 }
831
832 // Initialize the container to store the time offset.
833 // The number of parameters, 6, corresponds to the number of lines to be read
834 // from sTGC_Digitization_timeOffsetStrip.dat.
835 // Setting the default offset to 0 ns.
836
837 // Read the input file
838 std::string line;
839 size_t index{0};
840 double value{0.0};
841 while (std::getline(ifs, line)) {
842 std::string key;
843 std::istringstream iss(line);
844 iss >> key;
845 if (key.compare("strip") == 0) {
846 iss >> index >> value;
847 if (index >= m_timeOffsetStrip.size()) continue;
848 m_timeOffsetStrip.at(index) = value;
849 }
850 }
851 return StatusCode::SUCCESS;
852}

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level lvl)
inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

29{
30 m_lvl = lvl;
31}

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
mutableprivateinherited

Messaging initialized (initMessaging)

Definition at line 141 of file AthMessaging.h.

◆ m_applyAsBuiltBLines

bool sTgcDigitMaker::m_applyAsBuiltBLines {false}
private

◆ m_channelTypes

int sTgcDigitMaker::m_channelTypes {3}
private

define offsets and widths of time windows for signals from wiregangs and strips.

The offsets are defined as relative time diffference with respect to the time after TOF and cable length corrections. Bunch crossing time is specified.

Definition at line 147 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

147{3}; // 1 -> strips, 2 -> strips+wires, 3 -> strips/wires/pads

◆ m_chargeAngularFactor

double sTgcDigitMaker::m_chargeAngularFactor {4.0}
private

◆ m_clusterProfile

std::array<double, 2> sTgcDigitMaker::m_clusterProfile {0.573, 1.092}
staticconstexprprivate

Definition at line 161 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

161{0.573, 1.092};

◆ m_doPadSharing

bool sTgcDigitMaker::m_doPadSharing {false}
private

◆ m_doTimeOffsetStrip

bool sTgcDigitMaker::m_doTimeOffsetStrip {false}
private

◆ m_gammaParameter

std::vector<GammaParameter> sTgcDigitMaker::m_gammaParameter
private

◆ m_idHelperSvc

const Muon::IMuonIdHelperSvc* sTgcDigitMaker::m_idHelperSvc {nullptr}
private

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }
mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

135{ nullptr };

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }
mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

138{ MSG::NIL };

◆ m_meanGasGain

double sTgcDigitMaker::m_meanGasGain {5.e4}
private

Definition at line 149 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

149{5.e4}; // mean gain estimated from ATLAS note "ATL-MUON-PUB-2014-001"

◆ m_mostProbableArrivalTime

std::array<double, 5> sTgcDigitMaker::m_mostProbableArrivalTime {make_array<double, 5>(0.)}
private

Definition at line 133 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

constexpr std::array< T, N > make_array(const T &def_val)
Helper function to initialize in-place arrays with non-zero values.
Definition ArrayHelper.h:10

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls
mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels)

Definition at line 132 of file AthMessaging.h.

◆ m_nm

std::string AthMessaging::m_nm
privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

◆ m_posResAngular

double sTgcDigitMaker::m_posResAngular {0.305/m_StripResolution}
private

◆ m_posResIncident

double sTgcDigitMaker::m_posResIncident {1.}
private

◆ m_StripResolution

double sTgcDigitMaker::m_StripResolution {0.0949}
private

◆ m_theta

double sTgcDigitMaker::m_theta {10}
private

Definition at line 148 of file MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h.

148{10}; // theta=10 tuned to match cluster size and charge profile

◆ m_timeOffsetStrip

std::array<double, 6> sTgcDigitMaker::m_timeOffsetStrip {make_array<double, 6>(0.)}
private

The documentation for this class was generated from the following files: