Use linearity of the sin at leading order to check that the angular differences are either 0 or PI
230 {
232
234
235 static constexpr std::array<ChIndex ,4> statWithField{ChIndex::BIL, ChIndex::BML, ChIndex::BMS, ChIndex::BOL};
237 std::ranges::find(statWithField,
chIndex) != statWithField.end();
238
239
241
242
243 const TrkDriftCircleMath::Line&
line = segment.
line();
244
245 ATH_MSG_DEBUG(
"New segment: chi2 " << segment.
chi2() <<
" ndof " << segment.
ndof() <<
" line " <<
line.position().x() <<
","
246 <<
line.position().y() <<
" phi " <<
line.phi() <<
" associated clusters "
248
249
251 Amg::Vector3D lroaddir = sInfo.globalTrans.linear() * roaddir2;
252
253
254 double lxroad = 0.;
255
256 if (hasPhiMeasurements) {
257
258
259 double sphi = 0.;
260 double cphi = lroaddir.x();
261
262 if (isEndcap) {
263 sphi = lroaddir.y();
264 lxroad = lroadpos.x() + (-lroadpos.y() +
line.position().x()) * cphi / absmax(sphi, std::numeric_limits<double>::min());
265 } else {
266 sphi = lroaddir.z();
267 lxroad = lroadpos.x() + (-lroadpos.z() +
line.position().y()) * cphi / absmax(sphi, std::numeric_limits<double>::min());
268 }
269
270 double shortestTubeLen = 1e9;
271
272 for (
const TrkDriftCircleMath::DCOnTrack& driftCircle : segment.
dcs()) {
274
275 const MdtDriftCircleOnTrack* riodc{
driftCircle.rot()};
278 double tubelen = 0.5 * riodc->
prepRawData()->detectorElement()->getActiveTubeLength(lay, tube);
279 if (tubelen < shortestTubeLen) shortestTubeLen = tubelen;
280 }
281
282 if (std::abs(lxroad) > shortestTubeLen) {
283 ATH_MSG_DEBUG(
"coordinates far outside chamber! using global position of first hit ");
284 if (lxroad < 0.) shortestTubeLen *= -1.;
285 lxroad = shortestTubeLen;
286 }
287 } else {
288 lxroad = (sInfo.globalTrans * mdts[0]->prepRawData()->detectorElement()->surface(mdts[0]->
identify()).center()).x();
289 }
290
291
293
294
296
297
299 surfaceTransform.pretranslate(gpos);
300 double surfDim = 500.;
301 std::unique_ptr<Trk::PlaneSurface> surf = std::make_unique<Trk::PlaneSurface>(surfaceTransform, surfDim, surfDim);
302
303
305 double linephi =
line.phi();
306
307
308
310
311
312 std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>> > rioDistVec;
313
314
315 std::set<Identifier> deltaVec;
316 std::set<Identifier> outoftimeVec;
317
318 associateMDTsToSegment(gdir, segment, sInfo.geom, sInfo.globalTrans, sInfo.amdbTrans, deltaVec, outoftimeVec, rioDistVec, beta);
319 std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>>> garbage_collector;
320
321 TrkDriftCircleMath::DCSLHitSelector hitSelector;
322
324
325 TrkDriftCircleMath::DCSLFitter defaultFitter;
328 if (goodFit) {
330 std::abs(segment.
line().
x0() -
result.line().x0()) > 0.01 ||
331 std::abs(segment.
line().
y0() -
result.line().y0()) > 0.01) {
332
333 linephi =
result.line().phi();
334 lpos[1] =
result.line().position().x();
335 lpos[2] =
result.line().position().y();
336 gpos = sInfo.amdbTrans * lpos;
337
338
340 surfaceTransform.pretranslate(gpos);
341 surf = std::make_unique<Trk::PlaneSurface>(surfaceTransform, surfDim, surfDim);
342
343
345 }
346 }
347 }
348
349
350 Trk::LocalDirection segLocDir;
351 surf->globalToLocalDirection(gdir, segLocDir);
354 return nullptr;
355 }
356
357
361 if (std::min(std::abs(diff_phi), std::abs( std::abs(diff_phi) -
M_PI)) > 1.e-3 ||
362 std::min(std::abs(diff_prec), std::abs(std::abs(diff_prec) -
M_PI)) > 1.e-3) {
363 ATH_MSG_WARNING(
" ALARM updated angles wrong: diff phi " << diff_phi <<
" prec " << diff_prec <<
" phi rdir " << roaddir2.phi()
364 << " gdir " << gdir.phi() << " lphi " << linephi << " seg "
366 }
367
368
369 std::pair<std::pair<int, int>, bool> netaPhiHits =
371
372 if (rioDistVec.empty()){
374 return nullptr;
375
376 }
378 auto meas_for_fit = [&rioDistVec] () {
379 std::vector<const Trk::MeasurementBase*>
out{};
380 out.reserve(rioDistVec.size());
381 std::sort(rioDistVec.begin(), rioDistVec.end(), SortByDistanceToSegment());
382 for (
const std::pair<
double, std::unique_ptr<const Trk::MeasurementBase>>& ele : rioDistVec)
out.push_back(ele.second.get());
384 };
385
386
387 double dlocx{1000.}, dangleXZ{1000.}, qoverp{-99999.}, dqoverp{-99999.};
388 bool hasMeasuredCoordinate = false;
390 ATH_MSG_DEBUG(
" distance between first and last phi hit sufficient to perform 4D fit: phi " << gdir.phi() <<
" theta "
391 << gdir.theta());
392
394
395 if (track) {
396 if (isCurvedSegment &&
track->perigeeParameters() &&
track->perigeeParameters()->covariance()) {
399 }
400 hasMeasuredCoordinate = true;
401
403 updatedCov.setZero();
404 m_segmentFitter->updateSegmentParameters(*track, *surf, segLocPos, segLocDir, updatedCov);
408 } else {
410 }
411
413 surf->localToGlobal(segLocPos, gdir, gpos);
414 surf->localToGlobalDirection(segLocDir, gdir);
415
416 if (
track->measurementsOnTrack() && rioDistVec.size() !=
track->measurementsOnTrack()->size()) {
417 if (
track->measurementsOnTrack()->empty()) {
419 return nullptr;
420 }
421 ATH_MSG_DEBUG(
" ROT vector size changed after fit, updating ");
422 garbage_collector = std::move(rioDistVec);
423 rioDistVec.reserve(
track->measurementsOnTrack()->size());
425 for (
const Trk::TrackStateOnSurface* tsit : *
track->trackStateOnSurfaces()) {
427 if (!pars) continue;
428 if (!firstPars) firstPars =
pars;
429
430
431 const Trk::MeasurementBase* meas = tsit->measurementOnTrack();
432 if (!meas) continue;
435 rioDistVec.emplace_back(dist, meas->
uniqueClone());
436 }
437 }
438 } else {
440 netaPhiHits.second = false;
441 }
442 }
443
444
445
447 if (
updateSegmentPhi(gpos, gdir, segLocPos, segLocDir, *surf, meas_for_fit(), sInfo.phimin, sInfo.phimax)) {
448 surf->localToGlobal(segLocPos, gpos, gpos);
449 surf->localToGlobalDirection(segLocDir, gdir);
450 hasMeasuredCoordinate = true;
451 dlocx = 100.;
452 dangleXZ = 0.1;
453 }
454 }
455
456 if (msgLvl(MSG::DEBUG)) {
457 std::vector<const Trk::MeasurementBase*> debug_meas = meas_for_fit();
458 ATH_MSG_DEBUG(
" number of hits " << debug_meas.size() <<
" of which trigger " << netaPhiHits.first.first <<
" eta and "
459 << netaPhiHits.first.second << " phi ");
460 for (const Trk::MeasurementBase* mit : debug_meas) {
461 const Trk::RIO_OnTrack* rot = dynamic_cast<const Trk::RIO_OnTrack*>(mit);
462 if (rot) {
464 const MdtDriftCircleOnTrack* mdt = dynamic_cast<const MdtDriftCircleOnTrack*>(rot);
465 if (mdt)
467 <<
" radius " << std::setw(6) << mdt->
driftRadius() <<
" time " << std::setw(6) << mdt->
driftTime());
468 continue;
469 }
470 const CompetingMuonClustersOnTrack* crot = dynamic_cast<const CompetingMuonClustersOnTrack*>(mit);
471 if (crot) {
473 << " comp rot with hits " << crot->containedROTs().size());
474 continue;
475 }
477 }
478 }
479
480 std::vector<Identifier> holeVec =
calculateHoles(ctx, chid, gpos, gdir, hasMeasuredCoordinate, deltaVec, outoftimeVec, rioDistVec);
481
482
483 if (!outoftimeVec.empty()) holeVec.insert(holeVec.end(), std::make_move_iterator(outoftimeVec.begin()),
484 std::make_move_iterator(outoftimeVec.end()));
485 MuonSegmentQuality* quality =
new MuonSegmentQuality(segment.
chi2(), segment.
ndof(), holeVec);
486
487 const TrkDriftCircleMath::DCSLFitter* dcslFitter =
m_dcslFitProvider->getFitter();
492 } else {
494 }
495 }
496 bool hasFittedT0 = false;
497 double fittedT0{0}, errorFittedT0{1.};
499 hasFittedT0 = true;
502 errorFittedT0 = segment.
t0Error();
503 }
else if (dcslFitter &&
result.hasT0Shift()) {
504 fittedT0 =
result.t0Shift();
505 errorFittedT0 =
result.t0Error();
506 } else {
508 hasFittedT0 = false;
509 }
510 }
511
512 std::unique_ptr<MuonSegment> msegment;
513 if (isCurvedSegment) {
514 if (qoverp == -99999.) {
515 double charge = gpos.z() * std::tan(gdir.theta());
517
518 constexpr double BILALPHA(28.4366), BMLALPHA(62.8267), BMSALPHA(53.1259), BOLALPHA(29.7554);
521 dqoverp = M_SQRT2 * segment.
dtheta() / BILALPHA;
522 }
else if (
chIndex == ChIndex::BML) {
524 dqoverp = M_SQRT2 * segment.
dtheta() / BMLALPHA;
525 }
else if (
chIndex == ChIndex::BMS) {
527 dqoverp = M_SQRT2 * segment.
dtheta() / BMSALPHA;
528 }
else if (
chIndex == ChIndex::BOL) {
530 dqoverp = M_SQRT2 * segment.
dtheta() / BOLALPHA;
531 }
532 }
540
541 std::vector<Trk::DefinedParameter> defPars;
544 defPars.emplace_back(gdir.phi(),
Trk::phi);
545 defPars.emplace_back(gdir.theta(),
Trk::theta);
547 Trk::LocalParameters segLocPar(defPars);
548 msegment = std::make_unique<MuonSegment>(
549 std::move(segLocPar),
550 std::move(covMatrix),
551 surf.release(),
553 quality,
555 } else {
556
563 msegment =
564 std::make_unique<MuonSegment>(segLocPos,
565 segLocDir,
566 std::move(covMatrix),
567 surf.release(),
569 quality,
571 }
572
573 if (hasFittedT0) msegment->setT0Error(fittedT0, errorFittedT0);
574
575
577
578 if (msgLvl(MSG::DEBUG)) {
583 }
584 if (segmentQuality < 0) { return nullptr; }
585 return msegment;
586 }
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.14.0