133 {
134 const EventContext& ctx = Gaudi::Hive::currentContext();
135
136 if (measurements.front()->isVertex()) {
137 return;
138 }
139
143 }
144
145
147 double qOverP = fitParameters.qOverP();
149 if (p < 0.) {
152 }
153
154
155
156 bool energyGain = false;
157 bool haveDelimiter = false;
158 std::optional<TrackSurfaceIntersection>
intersection = std::nullopt;
159 int leadingScatterers = 0;
160 Trk::FitMeasurement* leadingScatterer = nullptr;
161 for (auto* measurement : measurements) {
162 if ((*measurement).isMaterialDelimiter()) {
163 haveDelimiter = true;
164 } else if ((*measurement).isScatterer()) {
165
166 if (!(*measurement).numberDoF()) {
167 ++leadingScatterers;
168 leadingScatterer = measurement;
169 } else {
170 if (std::abs(1. / (*measurement).qOverP()) > p)
171 energyGain = true;
172 break;
173 }
174 }
175 }
176
177
178 if (haveDelimiter && !leadingScatterers) {
179
180 haveDelimiter = false;
182 const Surface* firstMeasurementSurface = nullptr;
183 Trk::FitMeasurement* leadingOutlier = nullptr;
184 std::vector<Trk::FitMeasurement*> leadingOutliers;
185 const Surface* surface = nullptr;
186 for (auto* measurement : measurements) {
187 if ((*measurement).isMaterialDelimiter()) {
188 haveDelimiter = true;
189 endPosition = (*measurement).position();
190 surface = (*measurement).surface();
191 } else if ((*measurement).isPositionMeasurement()) {
192 if ((*measurement).isOutlier()) {
193 if (!firstMeasurementSurface)
194 leadingOutliers.push_back(measurement);
195 } else {
197 firstMeasurementSurface = (*measurement).surface();
200 }
201 if (!haveDelimiter)
202 continue;
203
204 }
205 } else if ((*measurement).isScatterer()) {
206 if (!surface)
207 continue;
208
209 if (std::abs(1. / (*measurement).qOverP()) > p)
210 energyGain = true;
211 break;
212 }
213 }
214
215
216
217
218
219 const Perigee perigee(fitParameters.position(), p * fitParameters.direction(),
220 charge, fitParameters.vertex());
221 bool haveMaterial = false;
222 const std::vector<const TrackStateOnSurface*>*
indetMaterial =
nullptr;
225
226 if (
msgLvl(MSG::VERBOSE)) {
230 " addLeadingMaterial: using extrapolateM from distance "
231 << direction.dot(fitParameters.position() - startPosition));
232 }
233
234
237 particleHypothesis, garbage);
238
239
241 std::vector<const TrackStateOnSurface*>::const_reverse_iterator
r =
244
245 if (!(**r).trackParameters() || !(**r).materialEffectsOnTrack() ||
247 (**r).trackParameters()->position() - endPosition) > 0.)
248 continue;
249
250 haveMaterial = true;
251 }
252 }
253 } else {
254 haveDelimiter = false;
255 }
256
257
258 if (haveDelimiter && !haveMaterial) {
259
261 " no leading material found with forward extrapolation"
262 << ", try again with back extrapolation ");
263
264
267
268 std::vector<const TrackStateOnSurface*>* indetMaterialF = nullptr;
269 const std::vector<const TrackStateOnSurface*>* indetMaterialR = nullptr;
272 const PlaneSurface plane(
intersection->position(), uvt);
278
282
283 if (indetMaterialR && !indetMaterialR->empty()) {
284 indetMaterialF = new std::vector<const TrackStateOnSurface*>;
285 indetMaterialF->reserve(indetMaterialR->size());
286
287 std::vector<const TrackStateOnSurface*>::const_reverse_iterator
r =
288 indetMaterialR->rbegin();
289 for (;
r != indetMaterialR->rend(); ++
r) {
290 indetMaterialF->push_back(*
r);
291 }
292
293 for (
r = indetMaterialF->rbegin();
r != indetMaterialF->rend(); ++
r) {
294
295 if (!(**r).trackParameters() || !(**r).materialEffectsOnTrack() ||
297 (**r).trackParameters()->position() - endPosition) > 0.)
298 continue;
299
300 haveMaterial = true;
301 }
303 indetMaterialF = nullptr;
304 }
305 }
306 delete indetMaterialR;
307 }
308
309
310
311 FitMeasurement* leadingMeas = nullptr;
313 std::vector<const TrackStateOnSurface*>::const_reverse_iterator
r =
316
317 if (!(**r).trackParameters() || !(**r).materialEffectsOnTrack() ||
318 intersection->direction().dot((**r).trackParameters()->position() -
319 endPosition) > 0.)
320 continue;
321
322
323 double eLoss = 0.;
324 const MaterialEffectsOnTrack* materialEffects =
325 dynamic_cast<const MaterialEffectsOnTrack*>(
326 (**r).materialEffectsOnTrack());
327 if (materialEffects) {
328 eLoss = std::abs(materialEffects->energyLoss()->deltaE());
329 if (energyGain)
330 eLoss = -eLoss;
331 }
332
333 if (leadingScatterers++ || !firstMeasurementSurface) {
335 std::optional<TrackSurfaceIntersection> const newIntersectionSTEP =
337 ctx, (**r).trackParameters()->associatedSurface(),
341 (**r).trackParameters()->associatedSurface(), *
intersection,
343 } else {
347 ctx, (**r).trackParameters()->associatedSurface(),
350 (**r).trackParameters()->associatedSurface(),
352 }
353
354
358 break;
359 }
360 leadingMeas =
361 new FitMeasurement((**r).materialEffectsOnTrack(),
364 } else {
365
366
367 for (std::vector<Trk::FitMeasurement*>::const_iterator l =
368 leadingOutliers.begin();
369 l != leadingOutliers.end(); ++l) {
370 leadingOutlier = leadingOutliers.back();
371 measurements.erase(
372 std::remove(measurements.begin(), measurements.end(), *l),
373 measurements.end());
374 }
375 leadingMeas = new FitMeasurement(
376 (**r).materialEffectsOnTrack()->thicknessInX0(), -eLoss,
379 firstMeasurementSurface);
380 leadingScatterer = leadingMeas;
381 }
383 leadingMeas->qOverP(
qOverP);
384
385
386
387 if (leadingOutlier) {
390 while (
391 leadingOutlier &&
393 radius) {
394 leadingOutliers.pop_back();
395 measurements.insert(measurements.begin(), leadingOutlier);
396 if (!leadingOutliers.empty()) {
397 leadingOutlier = leadingOutliers.back();
398 } else {
399 leadingOutlier = nullptr;
400 }
401 }
402 }
403
405
406 measurements.insert(measurements.begin(), leadingMeas);
407
408
409 if (materialEffects) {
412 } else {
414 }
415 }
416 }
417 }
418
419
420 if (leadingMeas) {
422 std::optional<TrackSurfaceIntersection> const newIntersectionSTEP =
428 } else {
434 :
m_intersector->intersectSurface(perigee.associatedSurface(),
436 }
437 } else {
439 }
442 }
443
444
445
447 leadingCovariance.setZero();
448 if (leadingScatterers) {
449 double leadingScattering = 0.;
450 double previousScattering = 0.;
451 double leadingX0Integral = 0.;
452 std::vector<Trk::FitMeasurement*>::reverse_iterator
m =
453 measurements.rbegin();
454 while (*m != leadingScatterer)
455 ++m;
456 for (;
m != measurements.rend(); ++
m) {
457 if (!(**m).isScatterer())
458 continue;
459 const MaterialEffectsOnTrack* materialEffects =
460 dynamic_cast<const MaterialEffectsOnTrack*>((**m).materialEffects());
461 if (!materialEffects)
462 continue;
463
464
465 leadingX0Integral += (**m).materialEffects()->thicknessInX0();
467 leadingScattering = leadingX0Integral * logTerm * logTerm;
468 const double scatteringAngle =
470 std::sqrt(leadingScattering - previousScattering);
471 previousScattering = leadingScattering;
472 (**m).scatteringAngle(scatteringAngle, leadingX0Integral);
473
474
475 double angleSquared = 1. / (**m).weight();
477 fitParameters.vertex())
478 .perp();
479 const double sinThetaSquared =
481 angleSquared *= angleSquared / sinThetaSquared;
482
483
484 leadingCovariance(0, 0) +=
deltaR *
deltaR * angleSquared;
485 leadingCovariance(0, 2) -=
deltaR * angleSquared;
486 leadingCovariance(2, 0) = leadingCovariance(0, 2);
487 leadingCovariance(2, 2) += angleSquared;
488
489
490 leadingCovariance(1, 1) +=
492 leadingCovariance(1, 3) +=
deltaR * angleSquared;
493 leadingCovariance(3, 1) = leadingCovariance(1, 3);
494 leadingCovariance(3, 3) += angleSquared * sinThetaSquared;
495 }
496 }
497
498
501 qOverP, leadingCovariance);
502 }
503
504 else {
505 fitParameters.update(fitParameters.position(), fitParameters.direction(),
506 qOverP, leadingCovariance);
507 }
508
509
511 if (!haveDelimiter)
514 }
515}
Scalar deltaR(const MatrixBase< Derived > &vec) const
#define ATH_MSG_VERBOSE(x)
double charge(const T &p)
#define AmgSymMatrix(dim)
bool msgLvl(const MSG::Level lvl) const
const TrackSurfaceIntersection & intersection(ExtrapolationType type) const
ToolHandle< IPropagator > m_stepPropagator
ToolHandle< IIntersector > m_intersector
Trk::MagneticFieldProperties m_stepField
void printMeasurements(std::vector< FitMeasurement * > &measurements) const
const std::vector< const TrackStateOnSurface * > * extrapolatedMaterial(const ToolHandle< IExtrapolator > &extrapolator, const TrackParameters ¶meters, const Surface &surface, PropDirection dir, const BoundaryCheck &boundsCheck, ParticleHypothesis particleHypothesis, Garbage_t &garbage) const
void indetMaterial(std::vector< FitMeasurement * > &measurements, ParticleHypothesis particleHypothesis, const TrackParameters &startParameters, Garbage_t &garbage) const
static void deleteMaterial(const std::vector< const TrackStateOnSurface * > *material, Garbage_t &garbage)
ToolHandle< IExtrapolator > m_extrapolator
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ FullField
Field is set to be realistic, but within a given Volume.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
DataModel_detail::iterator< DVL > remove(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end, const T &value)
Specialization of remove for DataVector/List.
1.14.0