140{
141
142
143 SG::ReadCondHandle<HTcalculator> readHandle{
m_HTReadKey,ctx};
144 const HTcalculator* HTcalc = (*readHandle);
145
146 if(HTcalc==nullptr) {
148 }
149
150
151 const InDet::TRTPIDNN* PIDNN = nullptr;
154 PIDNN = (*readHandlePIDNN);
155
156 if(PIDNN==nullptr) {
158 }
159 }
160
161
163
164
166 if (!perigee) return PIDvalues;
167
168
173
174
176 ATH_MSG_DEBUG(
" Track has negative theta or is VERY close to beampipe! "
177 "(tan(theta/2) < 0.0001). Returning default Pid values.");
178 return PIDvalues;
179 }
180
181 if (qOverP == 0.0) {
182 ATH_MSG_DEBUG (
" Track momentum infinite! (i.e. q/p = 0). Returning default Pid values.");
183 return PIDvalues;
184 }
185
186 double pTrk = fabs(1.0 / qOverP);
189
190
192
194 ATH_MSG_WARNING(
" Occupancy was outside allowed range! Returning default Pid values. Occupancy = "
196 return PIDvalues;
197 }
198
201 ATH_MSG_DEBUG (
"check---------------------------------------------------------------------------------------");
203 ATH_MSG_DEBUG (
"check---------------------------------------------------------------------------------------");
204
205
206 double pHTel_prod = 1.0;
207 double pHTpi_prod = 1.0;
208
209
210
211
212
213 std::vector<double> hit_HTMB;
214 std::vector<double> hit_gasType;
215 std::vector<double> hit_tot;
216 std::vector<double> hit_L;
217 std::vector<double> hit_rTrkWire;
218 std::vector<double> hit_HitZ;
219 std::vector<double> hit_HitR;
220 std::vector<double> hit_isPrec;
221
222 unsigned int nTRThits = 0;
223 unsigned int nTRThitsHTMB = 0;
224 unsigned int nXehits = 0;
225 unsigned int nArhits = 0;
226 unsigned int nPrecHits = 0;
227
228
229
231 if (not recoTrackStates) {
233
234 return PIDvalues;
235 }
236
239
240
241 for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
242
243 const Trk::MeasurementBase *measurement = (*tsosIter)->measurementOnTrack();
244 if (!measurement) continue;
245
246
247
248 const InDet::TRT_DriftCircleOnTrack* driftcircle = nullptr;
250 const Trk::RIO_OnTrack* tmpRio =
251 static_cast<const Trk::RIO_OnTrack*>(measurement);
253 driftcircle = static_cast<const InDet::TRT_DriftCircleOnTrack*>(tmpRio);
254 }
255 }
256
257 if (!driftcircle) continue;
258
259
260 bool isHTMB = (driftcircle->
prepRawData()->getWord() & 0x00020000) > 0;
261
262 nTRThits++;
263 if (isHTMB) nTRThitsHTMB++;
264 hit_HTMB.push_back(static_cast<double>(isHTMB));
265
266
267
268
269
270 Identifier DCid = driftcircle->
identify();
271
272
273 int TrtPart = 0;
274 if (abs(
m_trtId->barrel_ec(DCid)) == 2)
275 TrtPart = (
m_trtId->layer_or_wheel(DCid) < 6) ? 1 : 2;
276
277
278 int StrawLayer = 0;
279 if (TrtPart == 0) {
280
281 if (
m_trtId->layer_or_wheel(DCid) == 0) {
282 StrawLayer =
m_trtId->straw_layer(DCid);
283 }
else if (
m_trtId->layer_or_wheel(DCid) == 1) {
284 StrawLayer = 19 +
m_trtId->straw_layer(DCid);
285 } else {
286 StrawLayer = 19 + 24 +
m_trtId->straw_layer(DCid);
287 }
288 } else {
289
290 if (
m_trtId->layer_or_wheel(DCid) < 6) {
291 StrawLayer =
293 } else {
294 StrawLayer =
295 8 * (
m_trtId->layer_or_wheel(DCid) - 6) +
m_trtId->straw_layer(DCid);
296 }
297 }
298
299
300 double HitZ = 0.;
301 double HitR = 0.;
302 double rTrkWire = 0.;
303 bool hasTrackParameters = true;
304 if ((*tsosIter)->trackParameters()) {
305
307 HitR = gp.perp();
308 HitZ = gp.z();
309 rTrkWire = fabs((*tsosIter)->trackParameters()->parameters()[
Trk::driftRadius]);
310 if (rTrkWire > 2.2) rTrkWire = 2.175;
311 } else {
312
313 hasTrackParameters = false;
316 rTrkWire = 0;
317 }
318
319
321 hit_HitZ.push_back(HitZ);
322 hit_HitR.push_back(HitR);
323 hit_rTrkWire.push_back(rTrkWire);
326 }
327
328
329
330
331
332 int SL_max[3] = {73, 96, 64};
333 if (StrawLayer > SL_max[TrtPart] || StrawLayer < 0) {
334 ATH_MSG_WARNING(
" StrawLayer was outside allowed range! TrtPart = " << TrtPart <<
" SL = " << StrawLayer);
335 continue;
336 }
337
338 double ZRpos[3] = {fabs(HitZ), HitR, HitR};
339 double ZRpos_min[3] = { 0.0, 630.0, 630.0};
340 double ZRpos_max[3] = {720.0, 1030.0, 1030.0};
341 if (ZRpos[TrtPart] > ZRpos_max[TrtPart]) {
342 ATH_MSG_WARNING(
" ZRpos was above allowed range - adjusted! TrtPart = " << TrtPart <<
" ZRpos = " << ZRpos[TrtPart]);
343 ZRpos[TrtPart] = ZRpos_max[TrtPart] - 0.001;
344 }
345 if (ZRpos[TrtPart] < ZRpos_min[TrtPart]) {
346 ATH_MSG_WARNING(
" ZRpos was below allowed range - adjusted! TrtPart = " << TrtPart <<
" ZRpos = " << ZRpos[TrtPart]);
347 ZRpos[TrtPart] = ZRpos_min[TrtPart] + 0.001;
348 }
349
350
351
352
353
354
355
356 int GasType=0;
359 if ( stat==2 || stat==3 ) { GasType = 0; }
360 else if ( stat==1 || stat==4 ) { GasType = 1; }
361 else if ( stat==5 ) { GasType = 1; }
362 else if ( stat==6 ) { GasType = 1; }
363 else if ( stat==7 ) { GasType = 1;
364 }
365 else {
367 "getStatusHT = "
368 << stat
369 << ", must be 'Good(2)||Xenon(3)' or 'Dead(1)||Argon(4)' or "
370 "'Krypton(5)' or 'EmulatedArgon(6)' or 'EmulatedKr(7)'!");
371 throw std::exception();
372 }
373 }
374
376 << nTRThits << " TrtPart: " << TrtPart
377 << " GasType: " << GasType << " SL: " << StrawLayer
378 << " ZRpos: " << ZRpos[TrtPart] << " TWdist: " << rTrkWire
380
382
383 hit_gasType.push_back(static_cast<double>(GasType));
384 if (GasType == 0) {
385 nXehits++;
386 } else if (GasType == 1) {
387 nArhits++;
388 }
389
390
392 bool isPrec = false;
393 if (errDc < 1.0) {
394 isPrec = true;
395 nPrecHits++;
396 }
397 hit_isPrec.push_back(static_cast<double>(isPrec));
398 }
399
400
401
402
405 TrtPart,
406 GasType,
407 StrawLayer,
408 ZRpos[TrtPart],
409 rTrkWire,
411 hasTrackParameters);
414 TrtPart,
415 GasType,
416 StrawLayer,
417 ZRpos[TrtPart],
418 rTrkWire,
420 hasTrackParameters);
421
422 if (pHTel > 0.999 || pHTpi > 0.999 || pHTel < 0.001 || pHTpi < 0.001) {
424 << pHTel << " pHTpi = " << pHTpi
425 << " TrtPart: " << TrtPart << " SL: " << StrawLayer
426 << " ZRpos: " << ZRpos[TrtPart] << " TWdist: " << rTrkWire
428 continue;
429 }
430
431 if (pHTel > 0.80 || pHTpi > 0.50 || pHTel < 0.025 || pHTpi < 0.010) {
433 << pHTel << " pHTpi = " << pHTpi
434 << " TrtPart: " << TrtPart << " SL: " << StrawLayer
435 << " ZRpos: " << ZRpos[TrtPart] << " TWdist: " << rTrkWire
437 continue;
438 }
439
440
441 if (isHTMB) {pHTel_prod *= pHTel; pHTpi_prod *= pHTpi;}
442 else {pHTel_prod *= 1.0-pHTel; pHTpi_prod *= 1.0-pHTpi;}
443 ATH_MSG_DEBUG (
"check pHT(el): " << pHTel <<
" pHT(pi): " << pHTpi );
444
445 }
446
447
448
450
451
453
454 ATH_MSG_DEBUG (
"check nTRThits: " << nTRThits <<
" : " << nTRThitsHTMB
455 << " pHTel_prod: " << pHTel_prod
456 << " pHTpi_prod: " << pHTpi_prod
458
460 ctx,
461 &track,
462 true,
464
468 ctx,
469 &track,
470 false,
472
473 double dEdx_usedHits_noHTHits =
m_TRTdEdxTool->usedHits(ctx, &track,
false);
476
477
480
481
483 (limProbHT * limProbToT) /
484 ((limProbHT * limProbToT) + ((1.0 - limProbHT) * (1.0 - limProbToT)));
485
486
488
490 return PIDvalues;
491 }
492
493
494 std::map<std::string, std::map<std::string, double>> scalarInputs_NN = PIDNN->
getScalarInputs();
495 std::map<std::string, std::map<std::string, std::vector<double>>> vectorInputs_NN = PIDNN->
getVectorInputs();
496
497
498 double fAr = static_cast<double>(nArhits) / nTRThits;
499 double fHTMB = static_cast<double>(nTRThitsHTMB) / nTRThits;
500 double PHF = static_cast<double>(nPrecHits) / nTRThits;
501
502 if (!scalarInputs_NN.empty()) {
503 std::map<std::string, double>& trackVarMap = scalarInputs_NN.begin()->second;
511 storeNNVariable(trackVarMap,
"dEdx",
static_cast<double>(dEdx_noHTHits));
512 }
513
514 if (!vectorInputs_NN.empty()) {
515 std::map<std::string, std::vector<double>>& hitVarMap = vectorInputs_NN.begin()->second;
524 }
526
528 for (const auto& scalarInputs : scalarInputs_NN) {
530 for (const auto& variable : scalarInputs.second) {
532 }
533 }
534 for (const auto& vectorInputs : vectorInputs_NN) {
536 for (const auto& variable : vectorInputs.second) {
538 }
539 }
541
542 return PIDvalues;
543}
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_MSG_WARNING(x)
DataModel_detail::const_iterator< DataVector > const_iterator
float getProbHT(float pTrk, Trk::ParticleHypothesis hypothesis, int TrtPart, int GasType, int StrawLayer, float ZR, float rTrkAnode, float Occupancy, bool hasTrackPars) const
static float Limit(float prob)
const std::map< std::string, std::map< std::string, double > > & getScalarInputs() const
double evaluate(std::map< std::string, std::map< std::string, double > > &scalarInputs, std::map< std::string, std::map< std::string, std::vector< double > > > &vectorInputs) const
const std::map< std::string, std::map< std::string, std::vector< double > > > & getVectorInputs() const
double timeOverThreshold() const
returns time over threshold in ns for valid digits; zero otherwise
virtual const TRT_DriftCircle * prepRawData() const override final
returns the PrepRawData - is a TRT_DriftCircle in this scope
virtual const Amg::Vector3D & globalPosition() const override final
return the global position of this RIO_OnTrack
virtual const Trk::Surface & associatedSurface() const override final
returns the surface for the local to global transformation
static double calculateTrackLengthInStraw(const Trk::TrackStateOnSurface *trackState, const TRT_ID *identifier)
virtual bool type(MeasurementBaseType::Type type) const =0
Interface method checking the type.
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
Identifier identify() const
return the identifier -extends MeasurementBase
virtual bool rioType(RIO_OnTrackType::Type type) const =0
Method checking the Rio On Track type.
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Eigen::Matrix< double, 3, 1 > Vector3D
DataVector< const Trk::TrackStateOnSurface > TrackStates
@ eProbabilityBrem
Electron probability from Brem fitting (DNA).
@ TRTTrackOccupancy
TRT track occupancy.
@ eProbabilityComb
Electron probability from combining the below probabilities.
@ eProbabilityNN
Electron probability from NN.
@ eProbabilityNumberOfTRTHitsUsedFordEdx
Number of TRT hits used for dEdx measurement.
@ eProbabilityToT
Electron probability from Time-Over-Threshold (ToT) information.
@ TRTdEdx
dEdx from TRT ToT measurement.
@ eProbabilityHT
Electron probability from High Threshold (HT) information.
static const std::vector< float > eProbabilityDefault(numberOfeProbabilityTypes, 0.5)
ParametersBase< TrackParametersDim, Charged > TrackParameters