Use linearity of the sin at leading order to check that the angular differences are either 0 or PI
231 {
233
235
236 static constexpr std::array<ChIndex ,4> statWithField{ChIndex::BIL, ChIndex::BML, ChIndex::BMS, ChIndex::BOL};
238 std::ranges::find(statWithField,
chIndex) != statWithField.end();
239
240
242
243
244 const TrkDriftCircleMath::Line&
line = segment.
line();
245
246 ATH_MSG_DEBUG(
"New segment: chi2 " << segment.
chi2() <<
" ndof " << segment.
ndof() <<
" line " <<
line.position().x() <<
","
247 <<
line.position().y() <<
" phi " <<
line.phi() <<
" associated clusters "
249
250
252 Amg::Vector3D lroaddir = sInfo.globalTrans.linear() * roaddir2;
253
254
255 double lxroad = 0.;
256
257 if (hasPhiMeasurements) {
258
259
260 double sphi = 0.;
261 double cphi = lroaddir.x();
262
263 if (isEndcap) {
264 sphi = lroaddir.y();
265 lxroad = lroadpos.x() + (-lroadpos.y() +
line.position().x()) * cphi / absmax(sphi, std::numeric_limits<double>::min());
266 } else {
267 sphi = lroaddir.z();
268 lxroad = lroadpos.x() + (-lroadpos.z() +
line.position().y()) * cphi / absmax(sphi, std::numeric_limits<double>::min());
269 }
270
271 double shortestTubeLen = 1e9;
272
273 for (
const TrkDriftCircleMath::DCOnTrack& driftCircle : segment.
dcs()) {
275
276 const MdtDriftCircleOnTrack* riodc{
driftCircle.rot()};
279 double tubelen = 0.5 * riodc->
prepRawData()->detectorElement()->getActiveTubeLength(lay, tube);
280 if (tubelen < shortestTubeLen) shortestTubeLen = tubelen;
281 }
282
283 if (std::abs(lxroad) > shortestTubeLen) {
284 ATH_MSG_DEBUG(
"coordinates far outside chamber! using global position of first hit ");
285 if (lxroad < 0.) shortestTubeLen *= -1.;
286 lxroad = shortestTubeLen;
287 }
288 } else {
289 lxroad = (sInfo.globalTrans * mdts[0]->prepRawData()->detectorElement()->surface(mdts[0]->
identify()).center()).x();
290 }
291
292
294
295
297
298
300 surfaceTransform.pretranslate(gpos);
301 double surfDim = 500.;
302 std::unique_ptr<Trk::PlaneSurface> surf = std::make_unique<Trk::PlaneSurface>(surfaceTransform, surfDim, surfDim);
303
304
306 double linephi =
line.phi();
307
308
309
311
312
313 std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>> > rioDistVec;
314
315
316 std::set<Identifier> deltaVec;
317 std::set<Identifier> outoftimeVec;
318
319 associateMDTsToSegment(gdir, segment, sInfo.geom, sInfo.globalTrans, sInfo.amdbTrans, deltaVec, outoftimeVec, rioDistVec, beta);
320 std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>>> garbage_collector;
321
322 TrkDriftCircleMath::DCSLHitSelector hitSelector;
323
325
326 TrkDriftCircleMath::DCSLFitter defaultFitter;
329 if (goodFit) {
331 std::abs(segment.
line().
x0() -
result.line().x0()) > 0.01 ||
332 std::abs(segment.
line().
y0() -
result.line().y0()) > 0.01) {
333
334 linephi =
result.line().phi();
335 lpos[1] =
result.line().position().x();
336 lpos[2] =
result.line().position().y();
337 gpos = sInfo.amdbTrans * lpos;
338
339
341 surfaceTransform.pretranslate(gpos);
342 surf = std::make_unique<Trk::PlaneSurface>(surfaceTransform, surfDim, surfDim);
343
344
346 }
347 }
348 }
349
350
351 Trk::LocalDirection segLocDir;
352 surf->globalToLocalDirection(gdir, segLocDir);
355 return nullptr;
356 }
357
358
362 if (std::min(std::abs(diff_phi), std::abs( std::abs(diff_phi) -
M_PI)) > 1.e-3 ||
363 std::min(std::abs(diff_prec), std::abs(std::abs(diff_prec) -
M_PI)) > 1.e-3) {
364 ATH_MSG_WARNING(
" ALARM updated angles wrong: diff phi " << diff_phi <<
" prec " << diff_prec <<
" phi rdir " << roaddir2.phi()
365 << " gdir " << gdir.phi() << " lphi " << linephi << " seg "
367 }
368
369
370 std::pair<std::pair<int, int>, bool> netaPhiHits =
372
373 if (rioDistVec.empty()){
375 return nullptr;
376
377 }
379 auto meas_for_fit = [&rioDistVec] () {
380 std::vector<const Trk::MeasurementBase*>
out{};
381 out.reserve(rioDistVec.size());
382 std::sort(rioDistVec.begin(), rioDistVec.end(), SortByDistanceToSegment());
383 for (
const std::pair<
double, std::unique_ptr<const Trk::MeasurementBase>>& ele : rioDistVec)
out.push_back(ele.second.get());
385 };
386
387
388 double dlocx{1000.}, dangleXZ{1000.}, qoverp{-99999.}, dqoverp{-99999.};
389 bool hasMeasuredCoordinate = false;
391 ATH_MSG_DEBUG(
" distance between first and last phi hit sufficient to perform 4D fit: phi " << gdir.phi() <<
" theta "
392 << gdir.theta());
393
395
396 if (track) {
397 if (isCurvedSegment &&
track->perigeeParameters() &&
track->perigeeParameters()->covariance()) {
400 }
401 hasMeasuredCoordinate = true;
402
404 updatedCov.setZero();
405 m_segmentFitter->updateSegmentParameters(*track, *surf, segLocPos, segLocDir, updatedCov);
409 } else {
411 }
412
414 surf->localToGlobal(segLocPos, gdir, gpos);
415 surf->localToGlobalDirection(segLocDir, gdir);
416
417 if (
track->measurementsOnTrack() && rioDistVec.size() !=
track->measurementsOnTrack()->size()) {
418 if (
track->measurementsOnTrack()->empty()) {
420 return nullptr;
421 }
422 ATH_MSG_DEBUG(
" ROT vector size changed after fit, updating ");
423 garbage_collector = std::move(rioDistVec);
424 rioDistVec.reserve(
track->measurementsOnTrack()->size());
426 for (
const Trk::TrackStateOnSurface* tsit : *
track->trackStateOnSurfaces()) {
428 if (!pars) continue;
429 if (!firstPars) firstPars =
pars;
430
431
432 const Trk::MeasurementBase* meas = tsit->measurementOnTrack();
433 if (!meas) continue;
436 rioDistVec.emplace_back(dist, meas->
uniqueClone());
437 }
438 }
439 } else {
441 netaPhiHits.second = false;
442 }
443 }
444
445
446
448 if (
updateSegmentPhi(gpos, gdir, segLocPos, segLocDir, *surf, meas_for_fit(), sInfo.phimin, sInfo.phimax)) {
449 surf->localToGlobal(segLocPos, gpos, gpos);
450 surf->localToGlobalDirection(segLocDir, gdir);
451 hasMeasuredCoordinate = true;
452 dlocx = 100.;
453 dangleXZ = 0.1;
454 }
455 }
456
457 if (msgLvl(MSG::DEBUG)) {
458 std::vector<const Trk::MeasurementBase*> debug_meas = meas_for_fit();
459 ATH_MSG_DEBUG(
" number of hits " << debug_meas.size() <<
" of which trigger " << netaPhiHits.first.first <<
" eta and "
460 << netaPhiHits.first.second << " phi ");
461 for (const Trk::MeasurementBase* mit : debug_meas) {
462 const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(mit);
463 if (rot) {
465 const MdtDriftCircleOnTrack* mdt = dynamic_cast<const MdtDriftCircleOnTrack*>(rot);
466 if (mdt)
468 <<
" radius " << std::setw(6) << mdt->
driftRadius() <<
" time " << std::setw(6) << mdt->
driftTime());
469 continue;
470 }
471 const CompetingMuonClustersOnTrack* crot = dynamic_cast<const CompetingMuonClustersOnTrack*>(mit);
472 if (crot) {
474 << " comp rot with hits " << crot->containedROTs().size());
475 continue;
476 }
478 }
479 }
480
481 std::vector<Identifier> holeVec =
calculateHoles(ctx, chid, gpos, gdir, hasMeasuredCoordinate, deltaVec, outoftimeVec, rioDistVec);
482
483
484 if (!outoftimeVec.empty()) holeVec.insert(holeVec.end(), std::make_move_iterator(outoftimeVec.begin()),
485 std::make_move_iterator(outoftimeVec.end()));
486 MuonSegmentQuality* quality =
new MuonSegmentQuality(segment.
chi2(), segment.
ndof(), holeVec);
487
488 const TrkDriftCircleMath::DCSLFitter* dcslFitter =
m_dcslFitProvider->getFitter();
493 } else {
495 }
496 }
497 bool hasFittedT0 = false;
498 double fittedT0{0}, errorFittedT0{1.};
500 hasFittedT0 = true;
503 errorFittedT0 = segment.
t0Error();
504 }
else if (dcslFitter &&
result.hasT0Shift()) {
505 fittedT0 =
result.t0Shift();
506 errorFittedT0 =
result.t0Error();
507 } else {
509 hasFittedT0 = false;
510 }
511 }
512
513 std::unique_ptr<MuonSegment> msegment;
514 if (isCurvedSegment) {
515 if (qoverp == -99999.) {
516 double charge = gpos.z() * std::tan(gdir.theta());
518
519 constexpr double BILALPHA(28.4366), BMLALPHA(62.8267), BMSALPHA(53.1259), BOLALPHA(29.7554);
522 dqoverp = M_SQRT2 * segment.
dtheta() / BILALPHA;
523 }
else if (
chIndex == ChIndex::BML) {
525 dqoverp = M_SQRT2 * segment.
dtheta() / BMLALPHA;
526 }
else if (
chIndex == ChIndex::BMS) {
528 dqoverp = M_SQRT2 * segment.
dtheta() / BMSALPHA;
529 }
else if (
chIndex == ChIndex::BOL) {
531 dqoverp = M_SQRT2 * segment.
dtheta() / BOLALPHA;
532 }
533 }
541
542 std::vector<Trk::DefinedParameter> defPars;
545 defPars.emplace_back(gdir.phi(),
Trk::phi);
546 defPars.emplace_back(gdir.theta(),
Trk::theta);
548 Trk::LocalParameters segLocPar(defPars);
549 msegment = std::make_unique<MuonSegment>(
550 std::move(segLocPar),
551 std::move(covMatrix),
552 surf.release(),
554 quality,
556 } else {
557
564 msegment =
565 std::make_unique<MuonSegment>(segLocPos,
566 segLocDir,
567 std::move(covMatrix),
568 surf.release(),
570 quality,
572 }
573
574 if (hasFittedT0) msegment->setT0Error(fittedT0, errorFittedT0);
575
576
578
579 if (msgLvl(MSG::DEBUG)) {
584 }
585 if (segmentQuality < 0) { return nullptr; }
586 return msegment;
587 }
double charge(const T &p)
Gaudi::Property< bool > m_updatePhiUsingPhiHits
ToolHandle< IDCSLFitProvider > m_dcslFitProvider
Amg::Vector3D updateDirection(double linephi, const Trk::PlaneSurface &surf, const Amg::Vector3D &roaddir, bool isCurvedSegment) const
update the global direction, keeping the phi of the input road direction but using the local angle YZ
ToolHandle< IMuonSegmentSelectionTool > m_segmentSelectionTool
Gaudi::Property< bool > m_outputFittedT0
ToolHandle< IMuonSegmentFittingTool > m_segmentFitter
Gaudi::Property< bool > m_redo2DFit
Gaudi::Property< bool > m_refitParameters
bool updateSegmentPhi(const Amg::Vector3D &gpos, const Amg::Vector3D &gdir, Amg::Vector2D &segLocPos, Trk::LocalDirection &segLocDir, Trk::PlaneSurface &surf, const std::vector< const Trk::MeasurementBase * > &rots, double phimin, double phimax) const
void associateMDTsToSegment(const Amg::Vector3D &gdir, TrkDriftCircleMath::Segment &segment, const TrkDriftCircleMath::ChamberGeometry *multiGeo, const Amg::Transform3D &gToStation, const Amg::Transform3D &amdbToGlobal, std::set< Identifier > &deltaVec, std::set< Identifier > &outoftimeVec, std::vector< std::pair< double, std::unique_ptr< const Trk::MeasurementBase > > > &rioDistVec, double beta=1.) const
std::pair< std::pair< int, int >, bool > associateClustersToSegment(const TrkDriftCircleMath::Segment &segment, const Identifier &chid, const Amg::Transform3D &gToStation, ClusterVecPair &spVecs, double phimin, double phimax, std::vector< std::pair< double, std::unique_ptr< const Trk::MeasurementBase > > > &rioDistVec) const
static DataVector< const Trk::MeasurementBase > createROTVec(std::vector< std::pair< double, std::unique_ptr< const Trk::MeasurementBase > > > &rioDistVec)
std::vector< Identifier > calculateHoles(const EventContext &ctx, Identifier chid, const Amg::Vector3D &gpos, const Amg::Vector3D &gdir, bool hasMeasuredCoordinate, std::set< Identifier > &deltaVec, std::set< Identifier > &outoftimeVec, const std::vector< std::pair< double, std::unique_ptr< const Trk::MeasurementBase > > > &rioDistVec) const
virtual bool fit(Segment &result, const Line &line, const DCOnTrackVec &dcs, double t0Seed=-99999.) const
const HitSelection selectHitsOnTrack(const DCOnTrackVec &dcs) const
void hitsOnTrack(unsigned int hitsOnTrack)
unsigned int ndof() const
double angleXZ() const
access method for angle of local XZ projection
double angleYZ() const
access method for angle of local YZ projection
std::unique_ptr< MeasurementBase > uniqueClone() const
NVI Clone giving up unique pointer.
const Amg::Vector3D & momentum() const
Access method for the momentum.
const Amg::Vector3D & position() const
Access method for the position.
@ DCMathSegmentMakerCurved
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Affine3d Transform3D
ChIndex chIndex(const std::string &index)
convert ChIndex name string to enum
ChIndex
enum to classify the different chamber layers in the muon spectrometer
std::vector< DCOnTrack > DCOnTrackVec
@ loc2
generic first and second local coordinate
ParametersBase< TrackParametersDim, Charged > TrackParameters
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
const Identifier & identify(const UncalibratedMeasurement *meas)
Returns the associated identifier from the muon measurement.
1.16.1