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 (
auto 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 }
304 }
305 delete indetMaterialR;
306 }
307
308
309
310 FitMeasurement* leadingMeas = nullptr;
312 std::vector<const TrackStateOnSurface*>::const_reverse_iterator
r =
315
316 if (!(**r).trackParameters() || !(**r).materialEffectsOnTrack() ||
317 intersection->direction().dot((**r).trackParameters()->position() -
318 endPosition) > 0.)
319 continue;
320
321
322 double eLoss = 0.;
323 const MaterialEffectsOnTrack* materialEffects =
324 dynamic_cast<const MaterialEffectsOnTrack*>(
325 (**r).materialEffectsOnTrack());
326 if (materialEffects) {
327 eLoss = std::abs(materialEffects->energyLoss()->deltaE());
328 if (energyGain)
329 eLoss = -eLoss;
330 }
331
332 if (leadingScatterers++ || !firstMeasurementSurface) {
334 std::optional<TrackSurfaceIntersection> const newIntersectionSTEP =
336 ctx, (**r).trackParameters()->associatedSurface(),
340 (**r).trackParameters()->associatedSurface(), *
intersection,
342 } else {
346 ctx, (**r).trackParameters()->associatedSurface(),
349 (**r).trackParameters()->associatedSurface(),
351 }
352
353
357 break;
358 }
359 leadingMeas =
360 new FitMeasurement((**r).materialEffectsOnTrack(),
363 } else {
364
365
366 for (std::vector<Trk::FitMeasurement*>::const_iterator l =
367 leadingOutliers.begin();
368 l != leadingOutliers.end(); ++l) {
369 leadingOutlier = leadingOutliers.back();
370 measurements.erase(
371 std::remove(measurements.begin(), measurements.end(), *l),
372 measurements.end());
373 }
374 leadingMeas = new FitMeasurement(
375 (**r).materialEffectsOnTrack()->thicknessInX0(), -eLoss,
378 firstMeasurementSurface);
379 leadingScatterer = leadingMeas;
380 }
382 leadingMeas->qOverP(
qOverP);
383
384
385
386 if (leadingOutlier) {
389 while (
390 leadingOutlier &&
392 radius) {
393 leadingOutliers.pop_back();
394 measurements.insert(measurements.begin(), leadingOutlier);
395 if (!leadingOutliers.empty()) {
396 leadingOutlier = leadingOutliers.back();
397 } else {
398 leadingOutlier = nullptr;
399 }
400 }
401 }
402
404
405 measurements.insert(measurements.begin(), leadingMeas);
406
407
408 if (materialEffects) {
411 } else {
413 }
414 }
415 }
416 }
417
418
419 if (leadingMeas) {
421 std::optional<TrackSurfaceIntersection> const newIntersectionSTEP =
427 } else {
433 :
m_intersector->intersectSurface(perigee.associatedSurface(),
435 }
436 } else {
438 }
441 }
442
443
444
446 leadingCovariance.setZero();
447 if (leadingScatterers) {
448 double leadingScattering = 0.;
449 double previousScattering = 0.;
450 double leadingX0Integral = 0.;
451 std::vector<Trk::FitMeasurement*>::reverse_iterator
m =
452 measurements.rbegin();
453 while (m != measurements.rend() && *m != leadingScatterer) {
454 ++m;
455 }
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 += 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
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
for(size_t i=0;i< m_blockFillers.size();i++)
Fill one block.
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
const Amg::Vector3D & direction() const
Method to retrieve the direction at the Intersection.
@ FullField
Field is set to be realistic, but within a given Volume.
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
TrackSurfaceIntersection()=default
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.17.0