70 {
71
73 if(energyDeposit==0.) return {};
74
77 sTgcSimIdToOfflineId simToOffline{&idHelper};
78 const Identifier offlineId = simToOffline.
convert(hit.
sTGCId());
79 const Identifier layid = idHelper.
channelID(offlineId,
81 idHelper.
gasGap(offlineId),
83
87
88
90
91 const MuonGM::sTgcReadoutElement* detEl = condContainers.detMgr->getsTgcReadoutElement(layid);
92 if( !detEl ){ return {}; }
93
94
95
96
97
98 constexpr double length_tolerance = 0.01;
99
100
102 const Trk::PlaneSurface& SURF_WIRE = detEl->
surface(surfHash_wire);
103
104
106
108
110
111
114 Amg::Vector2D posOnSurf_wire{0.5 * (hitOnSurface_wire + pre_pos_wire_surf).block<2,1>(0,0)};
115
116
117
118
119
120
121
122
123
124
125
126
128
129
131 if ((wire_number < 1) || (wire_number > number_wires)) {
132
133
135 if ((wire_number < 1) || (wire_number > number_wires)) {
137 }
138 }
139
140
142 double dist_wire = ionization.distance;
143 if (dist_wire > 0.) {
144
145
146
148 if (ionization.posOnSegment.x() < ionization.posOnWire.x()) {
adjacent = -1;}
149
150
152
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;
157 ionization = std::move(ion_adj);
158 }
159 } else {
160 ATH_MSG_DEBUG(
"Failed to get the distance between the wire number = " << wire_number
162 << ". Number of wires = " << number_wires
164 return {};
165 }
166
167
168 posOnSurf_wire = ionization.posOnWire.block<2,1>(0,0);
169
171
173 <<
" posOnTrack: " <<
Amg::toString(ionization.posOnSegment, 3)
180
181
182
183
184
185
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.");
195 }
196
197
198 if (dist_wire > wire_pitch) {
199 return {};
200 }
201
202
204 const double par_kappa = gamParam.kParameter;
205 const double par_theta = gamParam.thetaParameter;
207
208 double gamma_mpv = (par_kappa - 1) * par_theta;
209
210 if (gamma_mpv < 0.) {gamma_mpv = 0.;}
211 const double t0_par = most_prob_time - gamma_mpv;
212
213
214
215
216
217
218 double digit_time = t0_par + CLHEP::RandGamma::shoot(condContainers.rndmEngine, par_kappa, 1/par_theta);
219
220
221
222
223
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
236 << ", time = " << digit_time
237 << ", k parameter = " << par_kappa
238 << ", theta parameter = " << par_theta
239 << ", most probable time = " << most_prob_time);
240
243 if (condContainers.efficiencies) {
244 const double efficiency = condContainers.efficiencies->getEfficiency(layid);
245
246 if (CLHEP::RandFlat::shoot(condContainers.rndmEngine,0.0,1.0) >
efficiency)
return {};
247 }
248
249
251
252
253 const double stripPropagationTime = 0.;
254
255 double sDigitTimeWire = digit_time;
256 double sDigitTimePad = sDigitTimeWire;
257 double sDigitTimeStrip = sDigitTimeWire + stripPropagationTime;
258
260
261
262
263
266
269
270
272 const Trk::PlaneSurface& SURF_STRIP = detEl->
surface(surfHash_strip);
273
276
278
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
293
294 int stripNumber = detEl->
stripNumber(posOnSurf_strip, newId);
295 if( stripNumber == -1 ){
296
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
302 if (stripNumber < 1) {
304 ATH_MSG_WARNING(
"Position on strip surface = (" << posOnSurf_strip.x() <<
", " << posOnSurf_strip.y() <<
")");
305 return {};
306 }
307 }
310 if(!
isValid && stripNumber != -1) {
312 return {};
313 }
314
316 const double stripHalfPitch = detEl->
channelPitch(newId)*0.5;
317
318
319
320
321
322
323
324 const double ionized_charge = (5.65E-6)*energyDeposit/CLHEP::keV;
325
326
327
328
329
330
331
332
333
335
336
337 const double total_charge =
gain*ionized_charge;
338
339
340
341 const double tan_theta = GLODIRE.perp()/GLODIRE.z();
342
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
349
350 const double norm = 0.5 * total_charge;
351
352
353
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.);
359 norm);
360
361
362 constexpr double tolerance_charge = 0.0005;
363
364 constexpr unsigned int max_neighbor = 10;
365
366
367 for (unsigned int iStrip = 0; iStrip <= max_neighbor; ++iStrip) {
368 int currentStrip = stripNumber + iStrip;
369 if (currentStrip > NumberOfStrips) break;
370
371
373 channelType, currentStrip,
isValid);
378 }
379
380
381
382 double x_relative = locpos.x() - peak_position;
383
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
386 if (
charge < tolerance_charge)
break;
387
388
389 double strip_time = sDigitTimeStrip;
390
391
393
394 int indexFromMiddleStrip = iStrip;
395 if ((iStrip > 0) && ((x_relative + (0.75 - iStrip * 2) * stripHalfPitch) < 0))
396 indexFromMiddleStrip = iStrip - 1;
397
399 }
400
402
404 << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
405 << ", neighbor index = " << iStrip
407 }
408 }
409
410
411 for (unsigned int iStrip = 1; iStrip <= max_neighbor; ++iStrip) {
412 int currentStrip = stripNumber - iStrip;
413 if (currentStrip < 1) break;
414
421 }
422
423
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
429 double strip_time = sDigitTimeStrip;
430
432 int indexFromMiddleStrip = ((x_relative + (iStrip * 2 - 0.75) * stripHalfPitch) > 0)? iStrip - 1 : iStrip;
434 }
435
437
439 << ", time = " << strip_time << ", time offset = " << strip_time-sDigitTimeStrip
440 << ", neighbor index = " << iStrip
442 }
443 }
444
445
448 return digits;
449 }
450
451
452
453
456
458 idHelper.
gasGap(layid), channelType, 1);
459
461 const Trk::PlaneSurface& SURF_PAD = detEl->
surface(surfHash_pad);
462
464 const Amg::Vector2D posOnSurf_pad(hitOnSurface_pad.x(), hitOnSurface_pad.y());
465
466
468
469 if(insideBounds) {
470 int padNumber = detEl->
stripNumber(posOnSurf_pad, PAD_ID);
471 if( padNumber == -1 ){
472
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;
478 padNumber = detEl->
stripNumber(new_position, PAD_ID);
479
480 if (padNumber < 1) {
482 ATH_MSG_WARNING(
"Position on pad surface = (" << posOnSurf_pad.x() <<
", " << posOnSurf_pad.y() <<
")");
483 return digits;
484 }
485 }
489
493 return digits;
494 }
495
499
500
501
502
503
504
505 double deltaX = halfPadWidthX - std::abs(
diff.x());
506 double deltaY= halfPadWidthY - std::abs(
diff.y());
509
510
511 if (deltaX < 0.)
deltaX = 0.1;
512 if (deltaY < 0.)
deltaY = 0.1;
513
516
517
518 int newPhi =
diff.x() > 0 ? padEtaPhi.second-1 : padEtaPhi.second+1;
519 int newEta =
diff.y() > 0 ? padEtaPhi.first+1 : padEtaPhi.first-1;
522
523 if (isNeighX && isNeighY && validEta && validPhi){
524
526 channelType, padEtaPhi.first, newPhi,
isValid);
528 channelType, newEta, padEtaPhi.second,
isValid);
530 channelType, newEta, newPhi,
isValid);
533
534
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
542 channelType, padEtaPhi.first, newPhi,
isValid);
543
544
545
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
553 channelType, newEta, padEtaPhi.second,
isValid);
554
555
556
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
565 addDigit(digits, newId, bctag, sDigitTimePad, 0.5*total_charge);
566 }
567
568 }
569 else if(padNumber != -1) {
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
579 return digits;
580 }
581
582
583
584
585
588
589
591 idHelper.
gasGap(layid), channelType, 1);
592
593
595
596 if(insideBounds) {
597
598
599 int wiregroupNumber = detEl->
stripNumber(posOnSurf_wire, WIREGP_ID);
600 if( wiregroupNumber == -1 ){
601
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
607 if (wiregroupNumber < 1) {
609 ATH_MSG_WARNING(
"Position on wire surface = (" << posOnSurf_wire.x() <<
", " << posOnSurf_wire.y() <<
")");
610 return digits;
611 }
612 }
613
614
616 channelType, wiregroupNumber,
isValid);
617
620 if(wiregroupNumber>=1&&wiregroupNumber<=NumberOfWiregroups) {
621 addDigit(digits, newId, bctag, sDigitTimeWire, total_charge);
622 }
623 }
624 else if (wiregroupNumber != -1){
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
632
633 return digits;
634}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
virtual const Trk::PlaneSurface & surface() const override
access to chamber surface (phi orientation), uses the first gas gap
double wirePitch(int gas_gap=1) const
single wire pitch.
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
double globalTime() const
const Amg::Vector3D & globalDirection() const
double kineticEnergy() const
const Amg::Vector3D & globalPosition() const
const HepMcParticleLink & particleLink() const
int particleEncoding() const
double depositEnergy() const
static double getPadChargeFraction(double distance)
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 m_chargeAngularFactor
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.
static constexpr std::array< double, 2 > m_clusterProfile
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)
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