 |
ATLAS Offline Software
|
#include <GlobalChi2Fitter.h>
|
| | GlobalChi2Fitter (const std::string &, const std::string &, const IInterface *) |
| |
| virtual | ~ GlobalChi2Fitter () |
| |
| virtual StatusCode | initialize () override |
| |
| virtual StatusCode | finalize () override |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const MeasurementSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const PrepRawDataSet &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const MeasurementSet &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
| |
| virtual Track * | alignmentFit (AlignmentCache &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=Trk::nonInteracting) const override |
| |
|
| Track * | fitIm (const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const |
| |
| Track * | myfit (const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const |
| |
| Track * | myfit_helper (Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const |
| |
| Track * | mainCombinationStrategy (const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const |
| |
| Track * | backupCombinationStrategy (const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const |
| |
| void | makeProtoState (Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const |
| |
| void | makeProtoStateFromMeasurement (Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const |
| |
| bool | processTrkVolume (Cache &, const Trk::TrackingVolume *tvol) const |
| |
| void | addMaterialUpdateTrajectory (Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * >> &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const |
| | Given layer information, probe those layers for scatterers and add them to a track. More...
|
| |
| void | addIDMaterialFast (const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const |
| | A faster strategy for adding scatter material to tracks, works only for inner detector tracks. More...
|
| |
| void | addMaterial (const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const |
| |
| std::vector< std::unique_ptr< TrackParameters > > | holesearchExtrapolation (const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const |
| | Helper method which performs an extrapolation with additional logic for hole search. More...
|
| |
| std::unique_ptr< const TrackParameters > | makePerigee (Cache &, const TrackParameters &, const ParticleHypothesis) const |
| |
| std::unique_ptr< const TrackParameters > | makeTrackFindPerigeeParameters (const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
| |
| std::unique_ptr< GXFTrackState > | makeTrackFindPerigee (const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
| |
| std::unique_ptr< Track > | makeTrack (const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
| |
| void | fillResidualsAndErrors (const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const |
| |
| void | tryToConverge (const Cache &cache, GXFTrajectory &trajectory, const int it) const |
| |
| void | updateSystemWithMaxBremPull (GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const |
| |
| void | fillDerivatives (GXFTrajectory &traj) const |
| |
| FitterStatusCode | runIteration (const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const |
| |
| FitterStatusCode | updateFitParameters (GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const |
| | Method to update peregee parameters, scattering angles, and brems. More...
|
| |
| void | updatePixelROTs (GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const |
| | Update the Pixel ROT using the current trajectory/local track parameters. More...
|
| |
| GXFTrajectory * | runTrackCleanerSilicon (const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const |
| |
| void | runTrackCleanerTRT (Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const |
| |
| PropagationResult | calculateTrackParametersPropagateHelper (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| | Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method. More...
|
| |
| PropagationResult | calculateTrackParametersPropagate (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| | Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possible holes. More...
|
| |
| std::vector< std::reference_wrapper< GXFTrackState > > | holeSearchStates (GXFTrajectory &trajectory) const |
| | Extracts a collection of track states which are important for hole search. More...
|
| |
| std::optional< GlobalChi2Fitter::TrackHoleCount > | holeSearchProcess (const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState >> &states) const |
| | Conduct a hole search between a list of states, possibly reusing existing information. More...
|
| |
| void | holeSearchHelper (const std::vector< std::unique_ptr< TrackParameters >> &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const |
| | Helper method for the hole search that does the actual counting of holes and dead modules. More...
|
| |
| FitterStatusCode | calculateTrackParameters (const EventContext &ctx, GXFTrajectory &, bool) const |
| |
| std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > | updateEnergyLoss (const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const |
| |
| void | calculateTrackErrors (GXFTrajectory &, Amg::SymMatrixX &, bool) const |
| |
| std::optional< TransportJacobian > | numericalDerivatives (const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const |
| |
| virtual int | iterationsOfLastFit () const |
| |
| virtual void | setMinIterations (int) |
| |
| bool | isMuonTrack (const Track &) const |
| |
| void | initFieldCache (const EventContext &ctx, Cache &cache) const |
| | Initialize a field cache inside a fit cache object. More...
|
| |
| void | throwFailedToGetTrackingGeomtry () const |
| |
| bool | ensureValidEntranceCalo (const EventContext &ctx, Cache &cache) const |
| |
| bool | ensureValidEntranceMuonSpectrometer (const EventContext &ctx, Cache &cache) const |
| |
| const TrackingGeometry * | trackingGeometry (Cache &cache, const EventContext &ctx) const |
| |
| const TrackingGeometry * | retrieveTrackingGeometry (const EventContext &ctx) const |
| |
|
| static std::optional< std::pair< Amg::Vector3D, double > > | addMaterialFindIntersectionDisc (Cache &cache, const DiscSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat) |
| | Find the intersection of a set of track parameters onto a disc surface. More...
|
| |
| static std::optional< std::pair< Amg::Vector3D, double > > | addMaterialFindIntersectionCyl (Cache &cache, const CylinderSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat) |
| | Find the intersection of a set of track parameters onto a cylindrical surface. More...
|
| |
| static void | addMaterialGetLayers (Cache &cache, std::vector< std::pair< const Layer *, const Layer * >> &layers, std::vector< std::pair< const Layer *, const Layer * >> &uplayers, const std::vector< std::unique_ptr< GXFTrackState >> &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat) |
| | Collect all possible layers that a given track could have passed through. More...
|
| |
| static void | makeTrackFillDerivativeMatrix (Cache &, GXFTrajectory &) |
| |
| static void | fillFirstLastMeasurement (Cache &cache, GXFTrajectory &trajectory) |
| |
| static void | fillBfromMeasurements (const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b) |
| |
| static void | fillAfromMeasurements (const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) |
| |
| static void | fillAfromScatterers (GXFTrajectory &trajectory, Amg::SymMatrixX &a) |
| |
| static bool | tryToWeightAfromMaterial (Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2) |
| |
| static void | compensatePhiWeights (Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) |
| |
| static void | calculateDerivatives (GXFTrajectory &) |
| |
| static bool | correctAngles (double &, double &) |
| |
|
| ToolHandle< IRIO_OnTrackCreator > | m_ROTcreator {this, "RotCreatorTool", "", ""} |
| |
| ToolHandle< IRIO_OnTrackCreator > | m_broadROTcreator {this, "BroadRotCreatorTool", "", ""} |
| |
| ToolHandle< IUpdator > | m_updator {this, "MeasurementUpdateTool", "", ""} |
| |
| ToolHandle< IExtrapolator > | m_extrapolator {this, "ExtrapolationTool", "Trk::Extrapolator/CosmicsExtrapolator", ""} |
| |
| ToolHandle< IMultipleScatteringUpdator > | m_scattool {this, "MultipleScatteringTool", "Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator", ""} |
| |
| ToolHandle< IEnergyLossUpdator > | m_elosstool {this, "EnergyLossTool", "Trk::EnergyLossUpdator/AtlasEnergyLossUpdator", ""} |
| |
| ToolHandle< IMaterialEffectsUpdator > | m_matupdator {this, "MaterialUpdateTool", "", ""} |
| |
| ToolHandle< IPropagator > | m_propagator {this, "PropagatorTool", "", ""} |
| |
| ToolHandle< INavigator > | m_navigator {this, "NavigatorTool", "Trk::Navigator/CosmicsNavigator", ""} |
| |
| ToolHandle< IResidualPullCalculator > | m_residualPullCalculator {this, "ResidualPullCalculatorTool", "Trk::ResidualPullCalculator/ResidualPullCalculator", ""} |
| |
| ToolHandle< Trk::ITrkMaterialProviderTool > | m_caloMaterialProvider {this, "CaloMaterialProvider", "Trk::TrkMaterialProviderTool/TrkMaterialProviderTool", ""} |
| |
| ToolHandle< IMaterialEffectsOnTrackProvider > | m_calotool {this, "MuidTool", "Rec::MuidMaterialEffectsOnTrackProvider/MuidMaterialEffectsOnTrackProvider", ""} |
| |
| ToolHandle< IMaterialEffectsOnTrackProvider > | m_calotoolparam {this, "MuidToolParam", "", ""} |
| |
| ToolHandle< IBoundaryCheckTool > | m_boundaryCheckTool {this, "BoundaryCheckTool", "", "Boundary checking tool for detector sensitivities" } |
| |
| SG::ReadCondHandleKey< TrackingGeometry > | m_trackingGeometryReadKey |
| |
| SG::ReadCondHandleKey< AtlasFieldCacheCondObj > | m_field_cache_key |
| |
| const AtlasDetectorID * | m_DetID = nullptr |
| |
| Gaudi::Property< bool > | m_signedradius {this, "SignedDriftRadius", true} |
| |
| Gaudi::Property< bool > | m_calomat {this, "MuidMat", false} |
| |
| Gaudi::Property< bool > | m_extmat {this, "ExtrapolatorMaterial", true} |
| |
| Gaudi::Property< bool > | m_fillderivmatrix {this, "FillDerivativeMatrix", false} |
| |
| Gaudi::Property< bool > | m_straightlineprop {this, "StraightLine", true} |
| |
| Gaudi::Property< bool > | m_extensioncuts {this, "TRTExtensionCuts", true} |
| |
| Gaudi::Property< bool > | m_sirecal {this, "RecalibrateSilicon", false} |
| |
| Gaudi::Property< bool > | m_trtrecal {this, "RecalibrateTRT", false} |
| |
| Gaudi::Property< bool > | m_kinkfinding {this, "KinkFinding", false} |
| |
| Gaudi::Property< bool > | m_decomposesegments {this, "DecomposeSegments", true} |
| |
| Gaudi::Property< bool > | m_getmaterialfromtrack {this, "GetMaterialFromTrack", true} |
| |
| Gaudi::Property< bool > | m_domeastrackpar {this, "MeasuredTrackParameters", true} |
| |
| Gaudi::Property< bool > | m_storemat {this, "StoreMaterialOnTrack", true} |
| |
| Gaudi::Property< bool > | m_redoderivs {this, "RecalculateDerivatives", false} |
| |
| Gaudi::Property< bool > | m_reintoutl {this, "ReintegrateOutliers", false} |
| |
| Gaudi::Property< bool > | m_acceleration {this, "Acceleration", false} |
| |
| Gaudi::Property< bool > | m_numderiv {this, "NumericalDerivs", false} |
| |
| Gaudi::Property< bool > | m_fiteloss {this, "FitEnergyLoss", false} |
| |
| Gaudi::Property< bool > | m_asymeloss {this, "AsymmetricEnergyLoss", true} |
| |
| Gaudi::Property< bool > | m_useCaloTG {this, "UseCaloTG", false} |
| |
| Gaudi::Property< bool > | m_rejectLargeNScat {this, "RejectLargeNScat", false} |
| |
| Gaudi::Property< bool > | m_createSummary {this, "CreateTrackSummary", true} |
| |
| Gaudi::Property< bool > | m_holeSearch {this, "DoHoleSearch", false} |
| |
| Gaudi::Property< double > | m_outlcut {this, "OutlierCut", 5.0} |
| |
| Gaudi::Property< double > | m_p {this, "Momentum", 0.0} |
| |
| Gaudi::Property< double > | m_chi2cut {this, "TrackChi2PerNDFCut", 1.e15} |
| |
| Gaudi::Property< double > | m_scalefactor {this, "TRTTubeHitCut", 2.5} |
| |
| Gaudi::Property< double > | m_minphfcut {this, "MinPHFCut", 0.} |
| |
| Gaudi::Property< int > | m_maxoutliers {this, "MaxOutliers", 10} |
| |
| Gaudi::Property< int > | m_maxit {this, "MaxIterations", 30} |
| |
| Gaudi::Property< int > | m_miniter {this, "MinimumIterations", 1} |
| |
| Gaudi::Property< int > | m_fixbrem {this, "FixBrem", -1} |
| |
| Gaudi::Property< int > | m_maxitPixelROT {this, "IterationsToRebuildPixelRots", 0} |
| |
| SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > | m_clusterSplitProbContainer {this, "ClusterSplitProbabilityName", "",""} |
| |
| Trk::Volume | m_idVolume |
| |
| std::array< std::atomic< unsigned int >, S_MAX_VALUE > m_fit_status | ATLAS_THREAD_SAFE = {} |
| |
Definition at line 156 of file GlobalChi2Fitter.h.
◆ FitterStatusType
| Enumerator |
|---|
| S_FITS | |
| S_SUCCESSFUL_FITS | |
| S_MAT_INV_FAIL | |
| S_NOT_ENOUGH_MEAS | |
| S_PROPAGATION_FAIL | |
| S_INVALID_ANGLES | |
| S_NOT_CONVERGENT | |
| S_HIGH_CHI2 | |
| S_LOW_MOMENTUM | |
| S_MAX_VALUE | |
Definition at line 176 of file GlobalChi2Fitter.h.
◆ GlobalChi2Fitter()
| Trk::GlobalChi2Fitter::GlobalChi2Fitter |
( |
const std::string & |
t, |
|
|
const std::string & |
n, |
|
|
const IInterface * |
p |
|
) |
| |
◆ ~ GlobalChi2Fitter()
◆ addIDMaterialFast()
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
For every track, we need to add its scatterers. That is to say, we need to determine which bits of non-active material the particle in question may have passed through and add them to the track. This is generally an expensive operation, but we can cut some corners if the track only consists of inner detector hits. Specifically, we can exploit the layer structure of the detector to find possible material hits more quickly and efficiently than using the standard material adding algorithm, which is addMaterial.
- Parameters
-
| [in,out] | cache | General cache object, as used everywhere. |
| [in,out] | trajectory | The current state of the track, respresented in the fitter's internal track representation. States may be added to this. |
| [in] | parameters | Starting parameters for the material addition step. |
| [in] | part | Standard representation of particle type, used to determine the behaviour of the particle as it traverses materials. |
Definition at line 3465 of file GlobalChi2Fitter.cxx.
3477 if (!caloEntranceIsValid) {
3486 cache.m_negdiscs.empty() &&
3487 cache.m_posdiscs.empty() &&
3488 cache.m_barrelcylinders.empty()
3503 cache.m_fastmat =
false;
3504 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3520 bool hasmat =
false;
3521 int indexoffset = 0;
3522 int lastmatindex = 0;
3523 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.trackStates();
3525 GXFTrackState *firstsistate =
nullptr;
3526 GXFTrackState *lastsistate =
nullptr;
3537 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3538 if (oldstates[
i]->materialEffects() !=
nullptr) {
3547 if (firstsistate ==
nullptr) {
3548 if (oldstates[
i]->trackParameters() ==
nullptr) {
3549 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3552 oldstates[
i]->associatedSurface(),
3555 trajectory.m_fieldprop,
3559 if (tmppar ==
nullptr)
return;
3561 oldstates[
i]->setTrackParameters(std::move(tmppar));
3563 firstsistate = oldstates[
i].get();
3565 lastsistate = oldstates[
i].get();
3573 if (lastsistate ==
nullptr) {
3574 throw std::logic_error(
"No track state");
3583 if (lastsistate->trackParameters() ==
nullptr) {
3584 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3587 lastsistate->associatedSurface(),
3589 trajectory.m_fieldprop,
3593 if (tmppar ==
nullptr)
return;
3595 lastsistate->setTrackParameters(std::move(tmppar));
3604 refpar = lastsistate->trackParameters();
3605 indexoffset = lastmatindex;
3607 refpar = firstsistate->trackParameters();
3620 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3621 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.upstreamMaterialLayers();
◆ addMaterial()
Definition at line 3634 of file GlobalChi2Fitter.cxx.
3641 if (refpar2 ==
nullptr) {
3644 const MeasurementBase *firstmuonhit =
nullptr;
3645 const MeasurementBase *lastmuonhit =
nullptr;
3646 const MeasurementBase *firstidhit =
3648 const MeasurementBase *lastidhit =
nullptr;
3649 const MeasurementBase *firsthit =
nullptr;
3650 const MeasurementBase *lasthit =
nullptr;
3651 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3652 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3653 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3654 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3655 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3656 bool matvec_used=
false;
3657 std::unique_ptr<TrackParameters> startmatpar1;
3658 std::unique_ptr<TrackParameters> startmatpar2;
3669 int npseudomuon1 = 0;
3670 int npseudomuon2 = 0;
3672 for (
auto & state :
states) {
3675 GXFMaterialEffects *meff = state->materialEffects();
3678 if (firstidhit ==
nullptr) {
3687 if (firsthit ==
nullptr) {
3688 firsthit = state->measurement();
3689 if (cache.m_acceleration) {
3690 if (
tp ==
nullptr) {
3694 state->associatedSurface(),
3700 if (
tp ==
nullptr) {
3704 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3710 lasthit = state->measurement();
3716 if (firstidhit ==
nullptr) {
3717 firstidhit = state->measurement();
3720 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3724 lastidhit = state->measurement();
3725 if (
tp !=
nullptr) {
3730 if (firstsiliconpar ==
nullptr) {
3731 firstsiliconpar =
tp;
3733 lastsiliconpar =
tp;
3745 if (firstmuonhit ==
nullptr) {
3746 firstmuonhit = state->measurement();
3747 if (
tp !=
nullptr) {
3751 lastmuonhit = state->measurement();
3752 if (
tp !=
nullptr) {
3758 if (meff->deltaE() == 0) {
3759 if (firstcalopar ==
nullptr) {
3760 firstcalopar = state->trackParameters();
3762 lastcalopar = state->trackParameters();
3764 if (firstmatpar ==
nullptr) {
3765 firstmatpar = state->trackParameters();
3770 std::unique_ptr<TrackParameters> refpar;
3773 if (trajectory.m_straightline &&
m_p != 0) {
3777 refpar = refpar2->associatedSurface().createUniqueTrackParameters(
3778 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3781 if (firstmatpar !=
nullptr) {
3786 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3790 const double mass = trajectory.mass();
3792 const AmgVector(5) & newpars = startmatpar2->parameters();
3793 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
3796 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3797 newpars[0], newpars[1], newpars[2], newpars[3],
3798 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3802 }
else if (trajectory.m_straightline &&
m_p != 0) {
3806 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3807 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3813 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3814 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3818 if ((firstidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3820 const DistanceSolution distsol = firstidhit->associatedSurface().straightLineDistanceEstimate(
3822 refpar->momentum().unit()
3827 if (
distance < 0 && distsol.numberOfSolutions() > 0 && !cache.m_acceleration) {
3828 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3830 const Surface *destsurf = &firstidhit->associatedSurface();
3831 std::unique_ptr<const TrackParameters> tmppar;
3833 if (firstmuonhit !=
nullptr) {
3835 if (caloEntranceIsValid) {
3838 *cache.m_caloEntrance,
3842 if (tmppar !=
nullptr) {
3843 destsurf = &tmppar->associatedSurface();
3848 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
3853 false, matEffects) );
3856 if (matvec && !matvec->empty()) {
3857 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3858 const MaterialEffectsBase *meb = (*matvec)[
i]->materialEffectsOnTrack();
3861 const MaterialEffectsOnTrack *meot =
static_cast < const MaterialEffectsOnTrack *
>(meb);
3862 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3863 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->clone() :
nullptr;
3864 meff->setSigmaDeltaE(0);
3865 matstates.push_back(std::make_unique<GXFTrackState>(
3867 std::unique_ptr<const TrackParameters>(newpars)
3877 if ((lastidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3878 const DistanceSolution distsol = lastidhit->associatedSurface().straightLineDistanceEstimate(
3880 refpar->momentum().unit()
3885 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
3886 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3887 const Surface *destsurf = &lastidhit->associatedSurface();
3888 std::unique_ptr<const TrackParameters> tmppar;
3889 std::unique_ptr<Surface> calosurf;
3890 if (firstmuonhit !=
nullptr) {
3892 if (caloEntranceIsValid) {
3895 *cache.m_caloEntrance,
3900 if (tmppar !=
nullptr) {
3901 const CylinderSurface *cylcalosurf =
nullptr;
3904 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3906 const DiscSurface *disccalosurf =
nullptr;
3909 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3911 if (cylcalosurf !=
nullptr) {
3913 const CylinderBounds & cylbounds = cylcalosurf->bounds();
3914 const double radius = cylbounds.r();
3915 const double hlength = cylbounds.halflengthZ();
3916 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
3917 }
else if (disccalosurf !=
nullptr) {
3918 const double newz = (
3919 disccalosurf->center().z() > 0 ?
3920 disccalosurf->center().z() - 1 :
3921 disccalosurf->center().z() + 1
3925 disccalosurf->center().x(),
3926 disccalosurf->center().y(),
3931 trans.translation() << newpos;
3933 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&disccalosurf->bounds());
3934 const double rmin = discbounds->rMin();
3935 const double rmax = discbounds->rMax();
3936 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3938 destsurf = calosurf.release();
3942 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
3944 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3945 matvec_used =
false;
3947 if (matvec && !matvec->empty()) {
3948 for (
const auto &
i : *matvec) {
3953 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
3954 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3955 if (cache.m_fiteloss && (meot->energyLoss() !=
nullptr)) {
3956 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
3959 if (matEffects ==
electron && cache.m_asymeloss) {
3960 meff->setDeltaE(-5);
3962 if (trajectory.numberOfTRTHits() == 0) {
3963 meff->setScatteringSigmas(0, 0);
3966 meff->setSigmaDeltaE(50);
3969 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
3971 matstates.push_back(std::make_unique<GXFTrackState>(
3973 std::unique_ptr<const TrackParameters>(newparams)
3985 if (cache.m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3988 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3992 firstmuonhit->associatedSurface(),
3997 if (calomeots.empty()) {
4002 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4003 if (lasthit == lastmuonhit) {
4004 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4007 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4010 calomeots[
i].associatedSurface(),
4013 trajectory.m_fieldprop,
4017 if (layerpar ==
nullptr) {
4022 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4025 lastcalopar = layerpar.get();
4029 const double qoverp = layerpar->parameters()[
Trk::qOverP];
4030 double qoverpbrem = 0;
4034 (firstmuonpar !=
nullptr) &&
4035 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4037 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4039 const double sign = (qoverp > 0) ? 1 : -1;
4040 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4043 const AmgVector(5) & newpar = layerpar->parameters();
4045 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4046 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4048 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4051 matstates.push_back(std::make_unique<GXFTrackState>(
4053 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4055 prevtrackpars = std::move(layerpar);
4060 firsthit == firstmuonhit &&
4061 (!cache.m_getmaterialfromtrack || lasthit == lastidhit)
4064 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4066 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4069 calomeots[
i].associatedSurface(),
4072 trajectory.m_fieldprop,
4076 if (layerpar ==
nullptr) {
4081 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4090 const double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4095 (lastmuonpar !=
nullptr) &&
4096 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4100 const double sign = (qoverpbrem > 0) ? 1 : -1;
4101 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4104 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4105 const AmgVector(5) & newpar = layerpar->parameters();
4107 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4108 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4112 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4117 if (lasthit == lastmuonhit && cache.m_extmat) {
4118 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4120 if (lastcalopar !=
nullptr) {
4122 if (msEntranceIsValid) {
4123 if (cache.m_msEntrance->inside(lastcalopar->position())) {
4126 *cache.m_msEntrance,
4130 if (muonpar1 !=
nullptr) {
4138 rot.col(2) = trackdir;
4140 trans.linear().matrix() << rot;
4141 trans.translation() << muonpar1->
position() - .1 * trackdir;
4142 PlaneSurface
const curvlinsurf(trans);
4144 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4152 if (curvlinpar !=
nullptr) {
4153 muonpar1 = std::move(curvlinpar);
4157 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->clone());
4161 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4164 DistanceSolution distsol;
4166 if (muonpar1 !=
nullptr) {
4167 distsol = lastmuonhit->associatedSurface().straightLineDistanceEstimate(
4175 if ((
distance > 0) and(distsol.numberOfSolutions() >
4176 0) and (firstmuonhit !=
nullptr)) {
4177 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4184 if (distsol.numberOfSolutions() == 1) {
4186 }
else if (distsol.numberOfSolutions() == 2) {
4188 std::abs(distsol.first()) < std::abs(distsol.second()) ?
4194 if (
distance < 0 && distsol.numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4195 if (firstmuonpar !=
nullptr) {
4198 if (trajectory.m_straightline &&
m_p != 0) {
4203 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4206 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4209 firstmuonhit->associatedSurface(),
4212 trajectory.m_fieldprop,
4216 if (tmppar !=
nullptr) {
4217 muonpar1 = std::move(tmppar);
4223 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4224 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4227 states.back()->associatedSurface(),
4231 matvec_used =
false;
4238 if (matvec && !matvec->empty()) {
4239 for (
int j = 0; j < (
int) matvec->size(); j++) {
4240 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4244 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4245 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4248 !trajectory.m_straightline &&
4249 (meot->energyLoss() !=
nullptr) &&
4250 std::abs(meff->deltaE()) > 25 &&
4251 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4253 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4256 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4258 matstates.push_back(std::make_unique<GXFTrackState>(
4260 std::unique_ptr<const TrackParameters>(newparams)
4270 if (firsthit == firstmuonhit && cache.m_extmat && (firstcalopar !=
nullptr)) {
4271 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4274 if (msEntranceIsValid) {
4275 if (cache.m_msEntrance->inside(firstcalopar->position())) {
4278 *cache.m_msEntrance,
4282 if (muonpar1 !=
nullptr) {
4290 rot.col(2) = trackdir;
4292 trans.linear().matrix() << rot;
4293 trans.translation() << muonpar1->
position() - .1 * trackdir;
4294 const PlaneSurface curvlinsurf(trans);
4296 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4304 if (curvlinpar !=
nullptr) {
4305 muonpar1 = std::move(curvlinpar);
4309 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->clone());
4313 DistanceSolution distsol;
4315 if (muonpar1 !=
nullptr) {
4316 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4324 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4326 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4327 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4330 states[0]->associatedSurface(),
4334 matvec_used =
false;
4336 if (matvec && !matvec->empty()) {
4337 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4339 for (
int j = 0; j < (
int) matvec->size(); j++) {
4340 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4342 if (meb !=
nullptr) {
4346 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4347 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4350 !trajectory.m_straightline &&
4351 (meot->energyLoss() !=
nullptr) &&
4352 std::abs(meff->deltaE()) > 25 &&
4353 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4355 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4359 (*matvec)[j]->trackParameters() !=
nullptr
4360 ? (*matvec)[j]->trackParameters()->clone()
4363 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4365 std::unique_ptr<const TrackParameters>(tmpparams)
4378 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4379 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4382 newstates.reserve(oldstates.size() + matstates.size());
4385 int firstlayerno = -1;
4387 if (cache.m_acceleration) {
4388 trajectory.addBasicState(std::move(oldstates[0]));
4394 for (
int i = cache.m_acceleration ? 1 : 0;
i < (
int) oldstates.size();
i++) {
4395 bool addlayer =
true;
4397 while (addlayer && layerno < (
int) matstates.size()) {
4399 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4401 const DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4402 layerpar->position(),
4403 layerpar->momentum().unit()
4408 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
4413 const double cylinderradius = layerpar->associatedSurface().bounds().r();
4414 const double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4416 if (trackimpact > cylinderradius - 5 *
mm) {
4422 if (
i == (
int) oldstates.size() - 1) {
4427 GXFMaterialEffects *meff = matstates[layerno]->materialEffects();
4429 if (meff->sigmaDeltaPhi() > .4 || (meff->sigmaDeltaPhi() == 0 && meff->sigmaDeltaE() <= 0)) {
4430 if (meff->sigmaDeltaPhi() > .4) {
4431 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4434 if (meff->sigmaDeltaPhi() == 0) {
4442 if (firstlayerno < 0) {
4443 firstlayerno = layerno;
4446 trajectory.addMaterialState(std::move(matstates[layerno]));
4448 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4449 const TrackingVolume *tvol =
m_navigator->volume(ctx,layerpar->position());
4450 const Layer *lay =
nullptr;
4452 if (tvol !=
nullptr) {
4453 lay = (tvol->closestMaterialLayer(layerpar->position(),layerpar->momentum().normalized())).
object;
4456 const MaterialProperties *matprop = lay !=
nullptr ? lay->fullUpdateMaterialProperties(*layerpar) :
nullptr;
4457 meff->setMaterialProperties(matprop);
4463 trajectory.addBasicState(std::move(oldstates[
i]));
4466 ATH_MSG_DEBUG(
"Total X0: " << trajectory.totalX0() <<
" total eloss: " << trajectory.totalEnergyLoss());
4468 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
◆ addMaterialFindIntersectionCyl()
Find the intersection of a set of track parameters onto a cylindrical surface.
See addMaterialFindIntersectionDisc for more information.
- Note
- This method can probably be replaced entirely by the straight line intersection method of the appropriate Surface subclass.
Definition at line 2995 of file GlobalChi2Fitter.cxx.
3009 const double *
pos = parforextrap.position().data();
3011 cache.m_field_cache.getFieldZR(
pos,
field);
3017 const double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3018 const double xc = parforextrap.position().x() -
r * sinphi;
3019 const double yc = parforextrap.position().y() +
r * cosphi;
3020 const double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
3021 const double z0 = parforextrap.position().z();
3022 const double d = xc * xc + yc * yc;
3023 const double rcyl = surf.bounds().r();
3024 double mysqrt = ((
r + rcyl) * (
r + rcyl) -
d) * (
d - (
r - rcyl) * (
r - rcyl));
3030 mysqrt = std::sqrt(mysqrt);
3031 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3032 double secondterm = (mysqrt * yc) / (2 *
d);
3033 const double x1 = firstterm + secondterm;
3034 const double x2 = firstterm - secondterm;
3035 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3036 secondterm = (mysqrt * xc) / (2 *
d);
3037 const double y1 = firstterm - secondterm;
3038 const double y2 = firstterm + secondterm;
3039 double x = parforextrap.position().x();
3040 double y = parforextrap.position().y();
3041 const double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3042 const double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3044 if (dist1 < dist2) {
3052 const double phi1 = std::atan2(
y - yc,
x - xc);
3055 const double delta_z =
r * deltaphi / tantheta;
3056 const double z =
z0 + delta_z;
3060 if (std::abs(
z - surf.center().z()) > surf.bounds().halflengthZ()) {
3069 const double costracksurf = std::abs(normal.unit().dot(trackdir));
3071 return std::make_pair(
intersect, costracksurf);
◆ addMaterialFindIntersectionDisc()
Find the intersection of a set of track parameters onto a disc surface.
Calculates the intersection from a point and momentum in space onto a disc surface which represents a disc-shaped layer in the detector. The position of the intersection can be used to find materials in that layer at that position.
- Parameters
-
| [in] | cache | The standard GX2F cache. |
| [in] | surface | The surface to intersect with. |
| [in] | param1 | The main track parameters to calculate the intersection from. |
| [in] | param2 | A secondary set of parameters used for electrons. The purpose of this is not known to us at this time. |
| [in] | mat | A particle hypothesis describing the behaviour of the particle. |
- Returns
- Nothing if the intersection failed (i.e. there was no intersection), otherwise both an intersection positition as well as the angle of inflection.
- Note
- This method can probably be replaced entirely by the straight line intersection method of the appropriate Surface subclass.
Definition at line 2950 of file GlobalChi2Fitter.cxx.
2962 const double *
pos = parforextrap.position().data();
2964 cache.m_field_cache.getFieldZR(
pos,
field);
2972 const double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
2973 const double xc = parforextrap.position().x() -
r * sinphi;
2974 const double yc = parforextrap.position().y() +
r * cosphi;
2975 const double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
2976 const double z0 = parforextrap.position().z();
2977 const double delta_s = (surf.center().z() -
z0) / costheta;
2978 const double delta_phi = delta_s * sintheta /
r;
2983 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&surf.bounds());
2985 if (
perp > discbounds->rMax() ||
perp < discbounds->rMin()) {
2989 const double costracksurf = std::abs(costheta);
2991 return std::make_pair(
intersect, costracksurf);
◆ addMaterialGetLayers()
Collect all possible layers that a given track could have passed through.
If we are to use layer information to determine possible scatterer hits, we must first gather those layers. That's what this method does. It looks for disc and barrel cylinder layers that the given track might have crossed and collects them into output vectors. One contains layers between states on the track, and the upstream layers lie before the first state of the track.
- Parameters
-
| [in,out] | cache | General cache object. |
| [out] | layers | Output vector for layers. |
| [out] | uplayers | Output vector for upstream layers, which lie before the first hit in the track. |
| [in] | states | A list of track states on the track. |
| [in] | first | The first track state. |
| [in] | last | The last track state. |
| [in] | refpar | Reference parameters from which to extrapolate. |
| [in] | hasmat | Are there any existing materials on this track? |
Definition at line 3285 of file GlobalChi2Fitter.cxx.
3298 upstreamlayers.reserve(5);
3305 const double firstz = firstsistate.trackParameters()->position().z();
3306 const double firstr = firstsistate.trackParameters()->position().perp();
3307 const double firstz2 = hasmat ? lastsistate.trackParameters()->position().z() : firstsistate.trackParameters()->position().z();
3308 const double firstr2 = hasmat ? lastsistate.trackParameters()->position().perp() : firstsistate.trackParameters()->position().perp();
3310 GXFTrackState *firststate = oldstates.front().get();
3311 GXFTrackState *laststate = oldstates.back().get();
3317 const double lastz = laststate->position().z();
3318 const double lastr = laststate->position().perp();
3320 const Layer *startlayer = firststate->associatedSurface().associatedLayer();
3321 const Layer *startlayer2 = hasmat ? lastsistate.associatedSurface().associatedLayer() :
nullptr;
3322 const Layer *endlayer = laststate->associatedSurface().associatedLayer();
3325 const double slope = (tantheta != 0) ? 1 / tantheta : 0;
3331 std::vector < const Layer *>::const_iterator
it;
3332 std::vector < const Layer *>::const_iterator itend;
3340 it = cache.m_posdiscs.begin();
3341 itend = cache.m_posdiscs.end();
3343 it = cache.m_negdiscs.begin();
3344 itend = cache.m_negdiscs.end();
3350 for (;
it != itend; ++
it) {
3355 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3363 const DiscBounds *discbounds =
static_cast<const DiscBounds *
> (&(*it)->surfaceRepresentation().bounds());
3368 if (discbounds->rMax() < firstr || discbounds->rMin() > lastr) {
3372 const double rintersect = firstr + ((*it)->surfaceRepresentation().center().z() - firstz) / slope;
3375 rintersect < discbounds->rMin() - 50 ||
3376 rintersect > discbounds->rMax() + 50
3386 if ((*
it) == endlayer) {
3398 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3401 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3408 (*
it) != startlayer &&
3409 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3410 (*it) == startlayer2)
3420 for (
const auto *barrelcylinder : cache.m_barrelcylinders) {
3424 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3431 const double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().r() - firstr) * slope;
3433 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3434 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3438 if (barrelcylinder == endlayer) {
3445 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3446 barrelcylinder == startlayer) {
3447 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3450 if (barrelcylinder != startlayer &&
3451 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3452 barrelcylinder == startlayer2)) {
3453 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
◆ addMaterialUpdateTrajectory()
Given layer information, probe those layers for scatterers and add them to a track.
This is the meat of the pudding, if you will. Given the information that we have about layers, go through them all and find any possible material hits that we need to add to the track.
- Parameters
-
| [in,out] | cache | General cache object. |
| [in,out] | track | The track object as it exists now in IR. |
| [in] | offset | The first state after any existing materials. |
| [in] | layers | The list of layers. |
| [in] | ref1 | The first set of reference parameters. |
| [in] | ref2 | The second set of reference parameters. |
| [in] | mat | The particle hypothesis describing the track behaviour. |
- Note
- Attentive readers may wonder why we pass this function a vector of layers, but not a vector of upstream layers. The reason for this is that the vector of upstream layers is also a member of the cache object.
Definition at line 3074 of file GlobalChi2Fitter.cxx.
3083 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3084 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(
states);
3095 for (
int i = 0;
i <= indexoffset;
i++) {
3096 trajectory.addBasicState(std::move(oldstates[
i]));
3105 for (
int i = indexoffset + 1;
i < (
int) oldstates.size();
i++) {
3106 const double rmeas = oldstates[
i]->position().perp();
3107 const double zmeas = oldstates[
i]->position().z();
3115 while (layerindex < (
int)
layers.size()) {
3117 double costracksurf = 0.0;
3132 const CylinderSurface *cylsurf =
static_cast<const CylinderSurface *
> (&
layer->surfaceRepresentation());
3138 if (oldstates[
i]->trackParameters() !=
nullptr) {
3139 const double rlayer = cylsurf->bounds().r();
3140 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->position().perp() - rlayer)) {
3141 parforextrap = oldstates[
i]->trackParameters();
3157 if (cylsurf->bounds().r() > rmeas)
break;
3165 const DiscSurface *discsurf =
static_cast<const DiscSurface *
> (&
layer->surfaceRepresentation());
3167 if (oldstates[
i]->trackParameters() !=
nullptr) {
3168 const double zlayer = discsurf->center().z();
3169 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->position().z() - zlayer)) {
3170 parforextrap = oldstates[
i]->trackParameters();
3181 if (std::abs(discsurf->center().z()) > std::abs(zmeas))
break;
3183 throw std::logic_error(
"Unhandled surface.");
3190 const MaterialProperties *matprop =
layer->layerMaterialProperties()->fullMaterial(
intersect);
3191 if (matprop ==
nullptr) {
3200 const double X0 = matprop->thicknessInX0();
3204 const double actualx0 =
X0 / costracksurf;
3205 const double de = -std::abs(
3206 (matprop->thickness() / costracksurf) *
3209 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3212 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3214 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3218 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3219 meff->setDeltaE(de);
3220 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3221 meff->setX0(actualx0);
3222 meff->setSurface(&
layer->surfaceRepresentation());
3223 meff->setMaterialProperties(matprop);
3229 std::unique_ptr<EnergyLoss> eloss;
3231 if (cache.m_fiteloss || (matEffects ==
electron && cache.m_asymeloss)) {
3232 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3234 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3239 if (eloss !=
nullptr) {
3240 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3244 if (matEffects ==
electron && cache.m_asymeloss) {
3245 meff->setDeltaE(-5);
3246 if (trajectory.numberOfTRTHits() == 0) {
3247 meff->setScatteringSigmas(0, 0);
3250 meff->setSigmaDeltaE(50);
3251 if (eloss !=
nullptr) {
3252 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3253 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3258 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3259 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3260 " sigma eloss: " << meff->sigmaDeltaE()
3267 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3269 std::unique_ptr<const TrackParameters>()
3272 trajectory.addMaterialState(std::move(matstate));
3281 trajectory.addBasicState(std::move(oldstates[
i]));
◆ alignmentFit()
Definition at line 1819 of file GlobalChi2Fitter.cxx.
1825 const EventContext& ctx = Gaudi::Hive::currentContext();
1829 delete alignCache.m_derivMatrix;
1830 alignCache.m_derivMatrix =
nullptr;
1832 delete alignCache.m_fullCovarianceMatrix;
1833 alignCache.m_fullCovarianceMatrix =
nullptr;
1834 alignCache.m_iterationsOfLastFit = 0;
1837 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1838 if(newTrack !=
nullptr){
1839 if(cache.m_derivmat.size() != 0)
1840 alignCache.m_derivMatrix =
new Amg::MatrixX(cache.m_derivmat);
1841 if(cache.m_fullcovmat.size() != 0)
1842 alignCache.m_fullCovarianceMatrix =
new Amg::MatrixX(cache.m_fullcovmat);
1843 alignCache.m_iterationsOfLastFit = cache.m_lastiter;
◆ backupCombinationStrategy()
Definition at line 1242 of file GlobalChi2Fitter.cxx.
1250 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1253 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1263 const auto *
const pParametersVector = indettrack->trackParameters();
1266 if (pParametersVector->size() > 1)
1267 firstidpar = (*pParametersVector)[1];
1269 firstidpar = pParametersVector->back();
1271 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1272 if ((firstidpar !=
nullptr) && (cache.m_caloEntrance !=
nullptr))
1276 if (lastidpar ==
nullptr) {
1277 lastidpar.reset(pParametersVector->back()->clone());
1280 std::unique_ptr < const TrackParameters > firstscatpar;
1281 std::unique_ptr < const TrackParameters > lastscatpar;
1282 std::unique_ptr < const TrackParameters > elosspar;
1287 lastidpar->position(),
1288 lastidpar->momentum(),
1290 lastidpar->position()
1297 calomeots[0].associatedSurface(),
1300 trajectory.m_fieldprop,
1303 if (!firstscatpar) {
1307 const std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1314 trajectory.m_fieldprop,
1324 const double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1333 const double newqoverp =
sign / std::sqrt(newp2);
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1345 calomeots[2].associatedSurface(),
1348 trajectory.m_fieldprop,
1358 calomeots[2].associatedSurface(),
1361 trajectory.m_fieldprop,
1372 calomeots[1].associatedSurface(),
1375 trajectory.m_fieldprop,
1383 const double sign = (elosspar->parameters()[
Trk::qOverP] < 0) ? -1 : 1;
1384 const double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1390 std::unique_ptr<const TrackParameters>
const tmppar(
1391 elosspar->associatedSurface().createUniqueTrackParameters(
1399 calomeots[0].associatedSurface(),
1402 trajectory.m_fieldprop,
1406 if (!firstscatpar) {
1411 for (; itStates != endState; ++itStates) {
1419 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1421 cache.m_idmat =
false;
1432 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1433 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1434 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1438 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1439 elossmeff->setSigmaDeltaE(sigmadp);
1442 elossmeff->setdelta_p(
dp);
1444 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1445 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1446 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1448 GXFTrackState *secondscatstate = trajectory.trackStates().back().get();
1449 const Surface *triggersurf1 =
nullptr;
1450 const Surface *triggersurf2 =
nullptr;
1454 bool seenmdt =
false;
1455 bool mdtbetweenphihits =
false;
1459 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1460 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1461 (!firstismuon ? ++itStates2 : --itStates2)
1464 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1469 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1470 const Surface *surf = &pMeasurement->associatedSurface();
1474 if (isCompetingRIOsOnTrack) {
1476 rot = &crot->rioOnTrack(0);
1479 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1482 if ((rot !=
nullptr) &&
m_DetID->
is_mdt(rot->identify()) && (triggersurf1 !=
nullptr)) {
1486 (rot !=
nullptr) && (
1492 const Amg::Vector3D measdir = surf->transform().rotation().col(0);
1493 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1494 const double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1496 const bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1500 (*itStates2)->trackParameters() !=
nullptr ?
1501 (*itStates2)->trackParameters()->position() :
1502 rot->globalPosition();
1503 if (triggersurf1 !=
nullptr) {
1504 triggerpos2 = thispos;
1505 triggersurf2 = surf;
1507 mdtbetweenphihits =
true;
1510 triggerpos1 = thispos;
1511 triggersurf1 = surf;
1517 double mdttrig1 = 999999;
1518 double mdttrig2 = 999999;
1519 const Surface *mdtsurf1 =
nullptr;
1520 const Surface *mdtsurf2 =
nullptr;
1523 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1524 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1525 (!firstismuon ? ++itStates2 : --itStates2)
1527 const Surface *surf =
nullptr;
1529 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1532 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1535 if (surf ==
nullptr) {
1538 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1542 if (isCompetingRioOnTrack) {
1544 rot = &crot->rioOnTrack(0);
1547 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1550 const bool thisismdt = rot and
m_DetID->
is_mdt(rot->identify());
1553 (*itStates2)->trackParameters() !=
nullptr ?
1554 (*itStates2)->trackParameters()->position() :
1555 pThisMeasurement->globalPosition();
1556 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).
mag() < mdttrig1) {
1557 mdttrig1 = (globpos - triggerpos1).
mag();
1560 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).
mag() < mdttrig2) {
1561 mdttrig2 = (globpos - triggerpos2).
mag();
1567 GXFTrackState * firstpseudostate =
nullptr;
1568 std::vector<GXFTrackState *> outlierstates;
1569 std::vector<GXFTrackState *> outlierstates2;
1571 outlierstates.reserve(10);
1572 outlierstates2.reserve(10);
1574 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1576 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1577 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1579 const bool isStraightLine =
1580 pMeasurement !=
nullptr ?
1587 (newpseudo ==
nullptr) && (
1588 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1589 std::abs(pMeasurement->globalPosition().z()) < 10000
1592 std::unique_ptr<const TrackParameters> par2;
1593 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1594 par2.reset((*itStates2)->trackParameters()->clone());
1598 *secondscatstate->trackParameters(),
1599 pMeasurement->associatedSurface(),
1602 trajectory.m_fieldprop,
1606 if (par2 ==
nullptr) {
1612 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1615 par2->associatedSurface()
1618 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1625 firstpseudo->setMeasurementErrors(
errors);
1626 firstpseudostate = firstpseudo.get();
1627 trajectory.addMeasurementState(std::move(firstpseudo));
1632 if (isPseudo && !firstismuon) {
1636 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1638 cache.m_idmat =
false;
1646 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1647 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1648 (mdtsurf1 !=
nullptr)
1650 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->transform());
1652 transf->translation() << triggerpos1;
1653 StraightLineSurface
const slsurf(*transf);
1657 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1661 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1668 pseudostate1->setMeasurementErrors(
errors);
1669 outlierstates2.push_back(pseudostate1.get());
1670 trajectory.addMeasurementState(std::move(pseudostate1));
1674 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1675 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1676 mdtbetweenphihits &&
1677 (mdtsurf2 !=
nullptr)
1679 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->transform());
1680 transf->translation() << triggerpos2;
1681 StraightLineSurface
const slsurf(*transf);
1685 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1689 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1696 pseudostate2->setMeasurementErrors(
errors);
1698 outlierstates2.push_back(pseudostate2.get());
1699 trajectory.addMeasurementState(std::move(pseudostate2));
1706 trajectory.trackStates().back()->measurementType() ==
TrackState::TGC ||
1708 trajectory.trackStates().back()->measurementType() ==
TrackState::RPC &&
1709 trajectory.trackStates().back()->measuresPhi()
1714 outlierstates.push_back(trajectory.trackStates().back().get());
1715 trajectory.setOutlier((
int) trajectory.trackStates().size() - 1,
true);
1720 trajectory.setNumberOfPerigeeParameters(0);
1724 trajectory.setPrefit(2);
1726 cache.m_matfilled =
true;
1727 const bool tmpacc = cache.m_acceleration;
1728 cache.m_acceleration =
false;
1730 const std::unique_ptr<Trk::Track> tmp_track(
1731 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1732 cache.m_acceleration = tmpacc;
1734 cache.m_matfilled =
false;
1737 trajectory.converged() &&
1738 std::abs(trajectory.residuals().tail<1>()(0) / trajectory.errors().tail<1>()(0)) > 10
1743 if (trajectory.converged()) {
1744 if (firstpseudostate !=
nullptr) {
1749 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1752 par2->associatedSurface()
1754 firstpseudostate->setMeasurement(std::move(newpseudo));
1755 firstpseudostate->setRecalibrated(
false);
1758 for (
int j = 0; j < (
int) trajectory.trackStates().size(); j++) {
1759 for (
const auto &
i : outlierstates2) {
1760 if (trajectory.trackStates()[j].get() ==
i) {
1761 trajectory.setOutlier(j,
true);
1765 for (
const auto &
i : outlierstates) {
1766 if (trajectory.trackStates()[j].get() ==
i) {
1767 trajectory.setOutlier(j,
false);
1773 itStates = (firstismuon ? beginStates2 : endState - 1);
1774 itStates != (firstismuon ? endState2 : beginStates - 1);
1775 (firstismuon ? ++itStates : --itStates)
1781 makeProtoState(cache, trajectory, *itStates, (firstismuon ? -1 : 0));
1785 trajectory.setPrefit(0);
1786 trajectory.setNumberOfPerigeeParameters(5);
1788 cache.m_matfilled =
false;
◆ calculateDerivatives()
| void Trk::GlobalChi2Fitter::calculateDerivatives |
( |
GXFTrajectory & |
trajectory | ) |
|
|
staticprivate |
Definition at line 8002 of file GlobalChi2Fitter.cxx.
8003 const int nstatesupstream = trajectory.numberOfUpstreamStates();
8004 const int nscatupstream = trajectory.numberOfUpstreamScatterers();
8005 const int nbremupstream = trajectory.numberOfUpstreamBrems();
8006 const int nscats = trajectory.numberOfScatterers();
8007 const int nperpars = trajectory.numberOfPerigeeParameters();
8008 const int nfitpars = trajectory.numberOfFitParameters();
8010 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8012 Matrix55 initialjac;
8013 initialjac.setZero();
8014 initialjac(4, 4) = 1;
8016 Matrix55 jacvertex(initialjac);
8018 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.numberOfScatterers(), initialjac);
8019 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.numberOfBrems(), initialjac);
8021 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
8022 GXFTrackState *prevstate =
nullptr;
8023 GXFTrackState *state =
nullptr;
8030 for (
const bool forward : {
false,
true}) {
8032 hit_begin = nstatesupstream;
8034 scatno = nscatupstream;
8035 bremno = nbremupstream;
8037 hit_begin = nstatesupstream - 1;
8039 scatno = trajectory.numberOfUpstreamScatterers() - 1;
8040 bremno = trajectory.numberOfUpstreamBrems() - 1;
8044 int hitno = hit_begin;
8045 forward ? (hitno < hit_end) : (hitno >= hit_end);
8046 hitno += (forward ? 1 : -1)
8049 state =
states[hitno].get();
8053 if (fillderivmat && state->derivatives().cols() != nfitpars) {
8054 state->derivatives().resize(5, nfitpars);
8055 state->derivatives().setZero();
8061 const int jmaxbrem = 4;
8063 if (hitno == (forward ? hit_end - 1 : 0)) {
8064 if (!fillderivmat) {
8072 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8074 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8075 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8076 jacvertex(4, 4) = jac(4, 4);
8086 jcnt = jmax - jmin + 1;
8088 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8091 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8093 i == scatno + (forward ? -1 : 1) &&
8094 prevstate !=
nullptr &&
8096 (!trajectory.prefit() || prevstate->materialEffects()->deltaE() == 0)
8098 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8099 jacscat[
i](4, 4) = jac(4, 4);
8101 calculateJac(jac, jacscat[
i], jmin, jmax);
8105 Eigen::MatrixXd & derivmat = state->derivatives();
8106 const int scatterPos = nperpars + 2 *
i;
8108 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
8114 jcnt = jmax - jmin + 1;
8116 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8119 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8121 i == bremno + (forward ? -1 : 1) &&
8123 prevstate->materialEffects() &&
8124 prevstate->materialEffects()->sigmaDeltaE() > 0
8126 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8127 jacbrem[
i](4, 4) = jac(4, 4);
8129 calculateJac(jac, jacbrem[
i], jmin, jmax);
8133 Eigen::MatrixXd & derivmat = state->derivatives();
8134 const int scatterPos = nperpars + 2 * nscats +
i;
8136 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
8140 calculateJac(jac, jacvertex, 0, 4);
8144 Eigen::MatrixXd & derivmat = state->derivatives();
8145 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8147 if (nperpars == 5) {
8148 derivmat.col(4).segment(0, 4) *= .001;
8149 derivmat(4, 4) = .001 * jacvertex(4, 4);
8155 (!trajectory.prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
8157 scatno += (forward ? 1 : -1);
8161 states[hitno]->materialEffects() &&
8162 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8164 bremno += (forward ? 1 : -1);
8167 prevstate =
states[hitno].get();
◆ calculateTrackErrors()
Definition at line 8174 of file GlobalChi2Fitter.cxx.
8182 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
8183 const int nstatesupstream = trajectory.numberOfUpstreamStates();
8185 GXFTrackState *prevstate =
nullptr;
8186 int i = nstatesupstream;
8187 for (
int j = 0; j < (
int)
states.size(); j++) {
8188 if (j < nstatesupstream) {
8195 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8196 if (stateno == 0 || stateno == nstatesupstream) {
8197 prevstate =
nullptr;
8200 std::unique_ptr<GXFTrackState> & state =
states[
index];
8201 if (state->materialEffects() !=
nullptr) {
8202 prevstate = state.get();
8206 if (!state->hasTrackCovariance()) {
8207 state->zeroTrackCovariance();
8209 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8211 if ((prevstate !=
nullptr) &&
8215 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8218 trackerrmat = jac * prevcov * jac.transpose();
8222 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8226 const MeasurementBase *measurement = state->measurement();
8227 const Amg::MatrixX & meascov = measurement->localCovariance();
8232 bool errorok =
true;
8233 for (
int i = 0;
i < 5;
i++) {
8236 && trackerrmat(
i,
i) > meascov(j, j)) {
8238 const double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8239 trackerrmat(
i,
i) = meascov(j, j);
8240 for (
int k = 0;
k < 5;
k++) {
8250 for (
int i = 0;
i < 5;
i++) {
8254 for (
int j = 0; j < 5; j++) {
8261 if (trajectory.m_straightline) {
8262 trackerrmat(4, 4) = 1
e-20;
8266 state->trackParameters();
8268 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8270 if (state->hasTrackCovariance()) {
8271 trkerrmat = (state->trackCovariance());
8273 trkerrmat = std::nullopt;
8276 const AmgVector(5) & tpars = tmptrackpar->parameters();
8277 std::unique_ptr<const TrackParameters> trackpar(
8278 tmptrackpar->associatedSurface().createUniqueTrackParameters(tpars[0],
8283 std::move(trkerrmat))
8285 state->setTrackParameters(std::move(trackpar));
8286 FitQualityOnSurface fitQual{};
8288 if (errorok && trajectory.nDOF() > 0) {
8289 fitQual =
m_updator->fullStateFitQuality(
8290 *state->trackParameters(),
8291 measurement->localParameters(),
8292 measurement->localCovariance()
8295 fitQual = FitQualityOnSurface(0, state->numberOfMeasuredParameters());
8298 state->setFitQuality(fitQual);
8300 prevstate = state.get();
◆ calculateTrackParameters()
Definition at line 7778 of file GlobalChi2Fitter.cxx.
7786 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7787 const int nstatesupstream = trajectory.numberOfUpstreamStates();
7788 const TrackParameters *prevtrackpar = trajectory.referenceParameters();
7789 std::unique_ptr<const TrackParameters> tmptrackpar;
7791 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7795 const DistanceSolution distsol = surf1.straightLineDistanceEstimate(
7796 prevtrackpar->position(), prevtrackpar->momentum().unit()
7803 distsol.numberOfSolutions() > 0 &&
7804 prevtrackpar != trajectory.referenceParameters()
7814 trajectory.m_fieldprop,
7821 (
rv.m_parameters !=
nullptr) &&
7822 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7828 if (
rv.m_parameters ==
nullptr) {
7829 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7830 " pos: " << prevtrackpar->position() <<
" destination surface: " << surf1);
7834 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7838 if (
rv.m_jacobian != std::nullopt) {
7840 states[hitno]->materialEffects() !=
nullptr &&
7841 states[hitno]->materialEffects()->deltaE() != 0 &&
7842 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7843 !trajectory.m_straightline
7845 const double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7846 const double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7847 const double mass = trajectory.mass();
7848 const double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7849 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7852 states[hitno]->setJacobian(*
rv.m_jacobian);
7853 }
else if (calcderiv) {
7858 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7860 if (meff !=
nullptr && hitno != 0) {
7861 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7862 surf, *meff, *
states[hitno]->trackParameters(), trajectory.mass(), -1
7865 if (std::holds_alternative<FitterStatusCode>(
r)) {
7866 return std::get<FitterStatusCode>(
r);
7869 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7870 prevtrackpar = tmptrackpar.get();
7872 prevtrackpar = currenttrackpar;
7876 prevtrackpar = trajectory.referenceParameters();
7878 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7881 const DistanceSolution distsol = surf.straightLineDistanceEstimate(prevtrackpar->position(), prevtrackpar->momentum().unit());
7885 if (
distance < 0 && distsol.numberOfSolutions() > 0 && prevtrackpar != trajectory.referenceParameters()) {
7894 trajectory.m_fieldprop,
7900 (
rv.m_parameters !=
nullptr) &&
7902 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7907 if (
rv.m_parameters ==
nullptr) {
7908 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7909 " pos: " << prevtrackpar->
7910 position() <<
" destination surface: " << surf);
7914 if (
rv.m_jacobian != std::nullopt) {
7916 states[hitno]->materialEffects() !=
nullptr &&
7917 states[hitno]->materialEffects()->deltaE() != 0 &&
7918 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7919 !trajectory.m_straightline
7921 const double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7922 const double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7923 const double mass = trajectory.mass();
7924 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7927 newp = std::sqrt(newp);
7930 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7933 states[hitno]->setJacobian(*
rv.m_jacobian);
7934 }
else if (calcderiv) {
7939 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7941 if (meff !=
nullptr) {
7942 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7943 surf, *meff, *
rv.m_parameters, trajectory.mass(), +1
7946 if (std::holds_alternative<FitterStatusCode>(
r)) {
7947 return std::get<FitterStatusCode>(
r);
7950 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7953 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7954 prevtrackpar =
states[hitno]->trackParameters();
◆ calculateTrackParametersPropagate()
Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possible holes.
This is a helper function for the calculateTrackParameters() method. Its purpose is to propagate from a set of track parameters onto a surface, finding the new track parameters at that surface and optionally the Jacobian and a list of possible holes.
This method uses another helper function, aptly called calculateTrackParametersPropagateHelper(), which wraps the actual propagator calls. What calculateTrackParametersPropagate() is call this method once and check the result. If the result is invalid, that is to say we didn't manage to extract a correct set of track parameters, we try again but with the propagation direction flipped.
If the calcderiv argument is set, this method will attempt to calculate the Jacobian as well as the new set of track parameters. This involves non-trivial logic which is abstracted away in the underlying helper function.
- Parameters
-
| [in] | ctx | An event context. |
| [in] | prev | The origin track parameters to start the propagation. |
| [in] | ts | The destination track state (in GX2F internal form). |
| [in] | propdir | The propagation direction. |
| [in] | bf | The magnetic field properties. |
| [in] | calcderiv | If set, calculate the derivative. |
| [in] | holesearch | If set, search for holes. |
- Returns
- An instance of PropagationResult, which is a struct with three members. Firstly, it contains a unique pointer to a set of track parameters, which are the track parameters at the destination track state following propagation. If these parameters are a nullpointer, that indicates a failure state. Secondly, if requested, the Jacobian is stored in this struct. This may be a nullptr if it was not requested or if the calculation of the Jacobian failed. Thirdly, it contains a vector of possible holes found between the start and end of the propagation. Since these hole states are not always necessary, they are wrapped in a std::optional type.
Definition at line 7752 of file GlobalChi2Fitter.cxx.
7761 PropagationResult
rv;
7764 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7767 if (
rv.m_parameters ==
nullptr) {
7768 propdir = invertPropdir(propdir);
7771 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
◆ calculateTrackParametersPropagateHelper()
Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method.
This method encapsulates some of the logic relating to passing or not passing a Jacobian matrix to make the calculateTrackParameters() a lot more readable.
For information about parameters see the IPropagator::propagateParameters() method, which this method almost directly wraps.
Definition at line 7713 of file GlobalChi2Fitter.cxx.
7722 std::unique_ptr<const TrackParameters>
rv;
7723 std::optional<TransportJacobian> jac{};
7734 if (
rv !=
nullptr && calcderiv) {
7739 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7745 return PropagationResult {
7748 std::move(extrapolation)
◆ compensatePhiWeights()
Definition at line 5953 of file GlobalChi2Fitter.cxx.
5958 const int nPerPars = trajectory.numberOfPerigeeParameters();
5959 std::size_t scatno = 0;
5961 for (
auto & state : trajectory.trackStates()) {
5962 const GXFMaterialEffects *meff = state->materialEffects();
5964 if (meff ==
nullptr || meff->sigmaDeltaPhi() == 0) {
5968 if (scatno >= cache.m_phiweight.size()) {
5970 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5971 throw std::range_error(
message.str());
5974 const bool isValidPlaneSurface =
5976 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5978 if (meff->deltaE() == 0 || isValidPlaneSurface) {
5979 const int scatNoIndex = 2 * scatno + nPerPars;
5980 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5981 cache.m_phiweight[scatno] = 1;
5988 if (meff->sigmaDeltaPhi() != 0) {
◆ correctAngles()
| bool Trk::GlobalChi2Fitter::correctAngles |
( |
double & |
phi, |
|
|
double & |
theta |
|
) |
| |
|
staticprivate |
◆ ensureValidEntranceCalo()
| bool Trk::GlobalChi2Fitter::ensureValidEntranceCalo |
( |
const EventContext & |
ctx, |
|
|
Cache & |
cache |
|
) |
| const |
|
private |
Definition at line 8534 of file GlobalChi2Fitter.cxx.
8535 if (cache.m_caloEntrance ==
nullptr) {
8539 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8547 if (cache.m_caloEntrance ==
nullptr) {
8552 return cache.m_caloEntrance !=
nullptr;
◆ ensureValidEntranceMuonSpectrometer()
| bool Trk::GlobalChi2Fitter::ensureValidEntranceMuonSpectrometer |
( |
const EventContext & |
ctx, |
|
|
Cache & |
cache |
|
) |
| const |
|
private |
Definition at line 8555 of file GlobalChi2Fitter.cxx.
8556 if (cache.m_msEntrance ==
nullptr) {
8560 cache.m_msEntrance =
geometry->trackingVolume(
"MuonSpectrometerEntrance");
8568 if (cache.m_msEntrance ==
nullptr) {
8573 return cache.m_msEntrance !=
nullptr;
◆ fillAfromMeasurements()
Definition at line 5750 of file GlobalChi2Fitter.cxx.
5755 const int nFitPars = trajectory.numberOfFitParameters();
5756 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5758 for (
int k = 0;
k < nFitPars;
k++) {
5759 for (
int l =
k;
l < nFitPars;
l++) {
5760 const int minMeas =
std::max(cache.m_firstmeasurement[
k], cache.m_firstmeasurement[
l]);
5761 const int maxMeas =
std::min(cache.m_lastmeasurement[
k], cache.m_lastmeasurement[
l]);
5764 for (
int measno = minMeas; measno < maxMeas; measno++) {
5765 a_kl += weightDeriv(measno,
k) * weightDeriv(measno,
l);
5768 a.fillSymmetric(
l,
k, a_kl);
◆ fillAfromScatterers()
Definition at line 5773 of file GlobalChi2Fitter.cxx.
5777 const int nFitPars = trajectory.numberOfFitParameters();
5778 const int nPerPars = trajectory.numberOfPerigeeParameters();
5779 const int nScatPars = 2 * trajectory.numberOfScatterers();
5780 const int nBrem = trajectory.numberOfBrems();
5781 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5784 const auto & scatSigmas = trajectory.scatteringSigmas();
5786 const int nMeas = (
int)
res.size();
5793 for (
int k = nPerPars;
k < nPerPars + nScatPars;
k += 2) {
5803 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5804 for (
int k = 4;
k < nFitPars;
k++) {
5806 k = nPerPars + nScatPars;
5809 for (
int l =
k;
l < nFitPars;
l++) {
5811 l = nPerPars + nScatPars;
5814 const double a_kl =
a(
l,
k) + weightDeriv(measno,
k) * weightDeriv(measno,
l);
5815 a.fillSymmetric(
l,
k, a_kl);
◆ fillBfromMeasurements()
Definition at line 5706 of file GlobalChi2Fitter.cxx.
5711 const int nFitPars = trajectory.numberOfFitParameters();
5712 const int nPerPars = trajectory.numberOfPerigeeParameters();
5713 const int nScatPars = 2 * trajectory.numberOfScatterers();
5714 const int nBrem = trajectory.numberOfBrems();
5715 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5720 const int nMeas = (
int)
res.size();
5722 for (
int k = 0;
k < nFitPars;
k++) {
5723 const int minMeasK = cache.m_firstmeasurement[
k];
5724 const int maxMeasK = cache.m_lastmeasurement[
k];
5731 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5732 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5742 if (
k == 4 ||
k >= nPerPars + nScatPars) {
5743 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5744 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
◆ fillDerivatives()
| void Trk::GlobalChi2Fitter::fillDerivatives |
( |
GXFTrajectory & |
traj | ) |
const |
|
private |
Definition at line 5497 of file GlobalChi2Fitter.cxx.
5502 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5506 const int nscatupstream = trajectory.numberOfUpstreamScatterers();
5507 const int nbremupstream = trajectory.numberOfUpstreamBrems();
5508 const int nscat = trajectory.numberOfScatterers();
5509 const int nbrem = trajectory.numberOfBrems();
5510 const int nperparams = trajectory.numberOfPerigeeParameters();
5512 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5515 const int nmeas = (
int) weightderiv.rows();
5517 for (std::unique_ptr<GXFTrackState> & state :
states) {
5520 const MeasurementBase *measbase = state->measurement();
5521 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5522 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5529 const double sinStereo =
5531 state->sinStereo() :
5533 const double cosStereo =
5535 std::sqrt(1 -
std::pow(sinStereo, 2)) :
5543 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5544 if (
i == 0 && sinStereo != 0) {
5545 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5547 return derivatives(
i, j);
5551 for (
int i = 0;
i < 5;
i++) {
5563 if (trajectory.numberOfPerigeeParameters() > 0) {
5564 const int cols = trajectory.m_straightline ? 4 : 5;
5566 if (
i == 0 && sinStereo != 0) {
5567 weightderiv.row(measno).head(
cols) =
5568 (derivatives.row(0).head(
cols) * cosStereo +
5569 sinStereo * derivatives.row(1).head(
cols)) /
5572 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5576 for (
int j = scatmin; j < scatmax; j++) {
5577 if (trajectory.prefit() == 1) {
5578 const int index = nperparams + j;
5581 const int index = nperparams + 2 * j;
5583 weightderiv(measno,
index + 1) = getThisDeriv(
i,
index + 1) /
error[measno];
5587 for (
int j = bremmin; j < bremmax; j++) {
5588 const int index = j + nperparams + 2 * nscat;
5595 double *
errors = state->measurementErrors();
5596 for (
int i = 0;
i < 5;
i++) {
5603 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5608 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5610 const double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5611 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5613 const double mass = .001 * trajectory.mass();
5615 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5619 auto multiplier = [] (
double mass,
double qOverP){
5622 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5623 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5625 if (trajectory.numberOfPerigeeParameters() > 0) {
5626 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5629 const auto bremNoBase = nperparams + 2 * nscat;
5630 if (bremno < nbremupstream) {
5631 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5632 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5633 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5636 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5637 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5638 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
◆ fillFirstLastMeasurement()
| void Trk::GlobalChi2Fitter::fillFirstLastMeasurement |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
trajectory |
|
) |
| |
|
staticprivate |
Definition at line 5646 of file GlobalChi2Fitter.cxx.
5650 const int nFitPars = trajectory.numberOfFitParameters();
5651 const int nPerPars = trajectory.numberOfPerigeeParameters();
5652 const int nScatPars = 2 * trajectory.numberOfScatterers();
5653 const int nBrem = trajectory.numberOfBrems();
5654 const int nUpstreamStates = trajectory.numberOfUpstreamStates();
5657 const int nMeas = (
int)
res.size();
5659 cache.m_firstmeasurement.resize(nFitPars);
5660 cache.m_lastmeasurement.resize(nFitPars);
5662 for (
int i = 0;
i < nPerPars;
i++) {
5663 cache.m_firstmeasurement[
i] = 0;
5664 cache.m_lastmeasurement[
i] = nMeas - nBrem;
5670 for (
int i = 0;
i < (
int) trajectory.trackStates().size();
i++) {
5671 const std::unique_ptr<GXFTrackState> & state = trajectory.trackStates()[
i];
5672 const GXFMaterialEffects *meff = state->materialEffects();
5674 if (meff ==
nullptr) {
5675 measno += state->numberOfMeasuredParameters();
5679 const int firstMeasurement =
i < nUpstreamStates ? 0 : measno;
5680 const int lastMeasurement =
i < nUpstreamStates ? measno : nMeas - nBrem;
5682 if (meff->sigmaDeltaTheta() != 0
5683 && (trajectory.prefit() == 0 || meff->deltaE() == 0)) {
5684 const int scatterPos = nPerPars + 2 * scatno;
5686 cache.m_firstmeasurement[scatterPos] = firstMeasurement;
5687 cache.m_lastmeasurement[scatterPos] = lastMeasurement;
5689 cache.m_firstmeasurement[scatterPos + 1] = firstMeasurement;
5690 cache.m_lastmeasurement[scatterPos + 1] = lastMeasurement;
5695 if (meff->sigmaDeltaE() > 0) {
5696 const int bremPos = nPerPars + nScatPars + bremno;
5698 cache.m_firstmeasurement[bremPos] = firstMeasurement;
5699 cache.m_lastmeasurement[bremPos] = lastMeasurement;
◆ fillResidualsAndErrors()
Definition at line 5073 of file GlobalChi2Fitter.cxx.
5084 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5104 const int nmeas = (
int)
res.size();
5105 const int nbrem = trajectory.numberOfBrems();
5106 const int nperpars = trajectory.numberOfPerigeeParameters();
5113 const int nidhits = trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits();
5114 const int nDOF = trajectory.nDOF();
5115 const bool doNewPseudoMeasurements = (
5119 std::abs((trajectory.prevchi2() - trajectory.chi2()) / nDOF) < 15 &&
5120 nidhits < trajectory.numberOfHits() &&
5121 (nperpars == 0 || nidhits > 0)
5132 double maxbrempull = -0.2;
5141 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5142 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5145 const MeasurementBase *measbase = state->measurement();
5156 doNewPseudoMeasurements &&
5158 !state->associatedSurface().isFree() &&
5159 !state->isRecalibrated()
5164 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5167 currenttrackpar->associatedSurface()
5170 state->setMeasurement(std::move(newpseudo));
5171 measbase = state->measurement();
5178 double *
errors = state->measurementErrors();
5180 for (
int i = 0;
i < 5;
i++) {
5200 res[measno] = residuals[
i];
5215 double *
errors = state->measurementErrors();
5216 for (
int i = 0;
i < 5;
i++) {
5229 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5231 const double deltaPhi = state->materialEffects()->deltaPhi();
5232 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5233 const double sigma2deltaPhi =
std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5234 const double deltaTheta = state->materialEffects()->deltaTheta();
5235 const double sigma2deltaTheta =
std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5237 if (trajectory.prefit() != 1) {
5238 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5239 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5241 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5246 deltaTheta * deltaTheta / sigma2deltaTheta
5255 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5256 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5258 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5259 const double pbrem = 1. / std::abs(qoverpbrem);
5260 const double p = 1. / std::abs(qoverp);
5261 const double mass = .001 * trajectory.mass();
5263 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5265 const double resMaterial = .001 * averagenergyloss -
energy + bremEnergy;
5266 res[nmeas - nbrem + bremno] = resMaterial;
5268 const double sigde = state->materialEffects()->sigmaDeltaE();
5269 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5270 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5272 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5273 error[nmeas - nbrem + bremno] = errorMaterial;
5284 if (state->materialEffects()->isKink()) {
5285 maxbrempull = -999999999;
5286 state_maxbrempull =
nullptr;
5290 cache.m_asymeloss &&
5292 trajectory.prefit() == 0 &&
5294 sigde != sigdepos &&
5297 const double elosspull = resMaterial / errorMaterial;
5299 if (trajectory.mass() > 100) {
5304 if (std::abs(elosspull) > 1) {
5305 if (elosspull < -1) {
5306 state->materialEffects()->setSigmaDeltaE(sigdepos);
5308 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5311 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5312 error[nmeas - nbrem + bremno] = errorMaterial;
5314 }
else if ((trajectory.numberOfTRTHits() == 0) ||
it >= 3) {
5324 !state->materialEffects()->isKink() && (
5325 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5329 bremno_maxbrempull = bremno;
5330 state_maxbrempull = state.get();
5331 maxbrempull = elosspull;
5340 trajectory.prefit() == 0 &&
5341 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5342 state->materialEffects()->isMeasuredEloss() &&
5343 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5348 std::vector<MaterialEffectsOnTrack> calomeots =
5353 parforcalo->associatedSurface(),
5361 if (calomeots.size() == 3) {
5362 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5363 const double newres = .001 * averagenergyloss -
energy + bremEnergy;
5364 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5366 const double oldPull = resMaterial / errorMaterial;
5367 const double newPull = newres / newerr;
5369 if (std::abs(newPull) < std::abs(oldPull)) {
5370 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5372 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5373 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5374 res[nmeas - nbrem + bremno] = newres;
5375 error[nmeas - nbrem + bremno] = newerr;
5379 state->materialEffects()->setMeasuredEloss(
false);
5389 for (
int imeas = 0; imeas < nmeas; imeas++) {
5390 if (
error[imeas] == 0) {
5400 trajectory.setPrevChi2(trajectory.chi2());
5401 trajectory.setChi2(
chi2);
◆ finalize()
| StatusCode Trk::GlobalChi2Fitter::finalize |
( |
| ) |
|
|
overridevirtual |
Definition at line 289 of file GlobalChi2Fitter.cxx.
292 if (m_fit_status[
S_FITS] > 0) {
295 <<
" track fits failed because of a matrix inversion failure");
297 <<
" tracks were rejected by the outlier logic");
299 <<
" track fits failed because of a propagation failure");
301 <<
" track fits failed because of an invalid angle (theta/phi)");
303 <<
" track fits failed because the fit did not converge");
305 <<
" tracks did not pass the chi^2 cut");
307 <<
" tracks were killed by the energy loss update");
310 return StatusCode::SUCCESS;
◆ fit() [1/6]
Definition at line 2401 of file GlobalChi2Fitter.cxx.
2408 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2413 GXFTrajectory trajectory;
2416 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2421 for (
const auto *itSet : rots) {
2422 if (itSet ==
nullptr) {
2423 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2429 std::unique_ptr<const TrackParameters> startpar(param.clone());
2432 matEffects ==
muon &&
2433 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0
2435 cache.m_matfilled =
true;
2436 trajectory.setPrefit(2);
2438 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2440 cache.m_matfilled =
false;
2442 if (!trajectory.converged()) {
2446 trajectory.setConverged(
false);
2447 const TrackParameters *firstpar = trajectory.trackStates()[0]->trackParameters();
2448 const TrackParameters *lastpar = trajectory.trackStates().back()->trackParameters();
2450 PerigeeSurface
const persurf(firstpar->position() - 10 * firstpar->momentum().unit());
2456 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2459 firstpar->associatedSurface()
2462 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2469 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2472 lastpar->associatedSurface()
2475 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2478 if (!trajectory.m_straightline) {
2479 trajectory.setPrefit(3);
2480 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2481 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2482 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2487 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2489 cache.m_matfilled =
true;
2491 if (!trajectory.converged()) {
2496 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2497 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2498 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2502 trajectory.setPrefit(0);
2505 firstpar = trajectory.trackStates().front()->trackParameters();
2510 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2513 firstpar->associatedSurface()
2516 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2520 trajectory.trackStates().front()->setMeasurementErrors(
errors);
2524 lastpar = trajectory.trackStates().back()->trackParameters();
2529 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2532 lastpar->associatedSurface()
2535 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2539 trajectory.trackStates().back()->setMeasurementErrors(
errors);
2543 std::unique_ptr<Track>
track;
2545 if (startpar !=
nullptr) {
2546 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2549 if (
track !=
nullptr) {
2552 cache.m_matfilled =
false;
◆ fit() [2/6]
Definition at line 2180 of file GlobalChi2Fitter.cxx.
2187 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2190 for (
const auto *prd : prds) {
2191 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2193 const PlaneSurface *plsurf =
nullptr;
2196 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2198 const StraightLineSurface *slsurf =
nullptr;
2201 slsurf =
static_cast < const StraightLineSurface *
>(&prdsurf);
2203 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2204 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2209 }
else if (slsurf !=
nullptr) {
2218 }
else if (plsurf !=
nullptr) {
2219 if (param.covariance() !=
nullptr) {
2241 if (rot !=
nullptr) {
2242 rots.push_back(rot);
2246 std::unique_ptr<Track>
track =
2247 fit(ctx, rots, param, runOutlier, matEffects);
2249 for (
const auto *rot : rots) {
◆ fit() [3/6]
Definition at line 2257 of file GlobalChi2Fitter.cxx.
2264 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2269 GXFTrajectory trajectory;
2272 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2279 if (minpar ==
nullptr) {
2280 minpar = *(inputTrack.trackParameters()->begin());
2289 const bool old_reintoutl = cache.m_reintoutl;
2290 cache.m_reintoutl =
false;
2291 const bool tmpasymeloss = cache.m_asymeloss;
2294 cache.m_asymeloss =
true;
2297 for (; itStates != endState; ++itStates) {
2301 !trajectory.trackStates().empty() &&
2302 (trajectory.trackStates().back()->materialEffects() !=
nullptr) &&
2303 trajectory.trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2305 trajectory.trackStates().back()->materialEffects()->setKink(
true);
2309 cache.m_reintoutl = old_reintoutl;
2311 for (
const auto & measBase : addMeasColl) {
2312 if (measBase ==
nullptr) {
2313 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2321 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2322 cache.m_asymeloss = tmpasymeloss;
2324 if (
track !=
nullptr) {
2325 const double oldqual =
2326 inputTrack.fitQuality()->numberDoF() != 0 ?
2327 inputTrack.fitQuality()->chiSquared() / inputTrack.fitQuality()->numberDoF() :
2330 const double newqual =
2331 track->fitQuality()->numberDoF() != 0 ?
2332 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() :
2335 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2339 track.reset(
nullptr);
2343 if (
track !=
nullptr) {
◆ fit() [4/6]
Definition at line 2353 of file GlobalChi2Fitter.cxx.
2360 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2364 for (
const auto *prd : prds) {
2365 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2368 std::unique_ptr<const TrackParameters>
const trackparForCorrect(
2369 hitparam->associatedSurface().createUniqueTrackParameters(
2384 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2387 if (rot !=
nullptr) {
2388 rots.push_back(rot);
2392 std::unique_ptr<Track>
track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2394 for (
const auto *rot : rots) {
◆ fit() [5/6]
Definition at line 1795 of file GlobalChi2Fitter.cxx.
1801 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1806 GXFTrajectory trajectory;
1809 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1814 return std::unique_ptr<Track>(
1815 fitIm(ctx, cache, inputTrack, runOutlier, matEffects));
◆ fit() [6/6]
Definition at line 316 of file GlobalChi2Fitter.cxx.
323 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
328 GXFTrajectory trajectory;
330 trajectory.m_straightline = (
331 !cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn()
338 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
339 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
341 bool measphi =
false;
343 for (
const auto *
i : *(muontrack->measurementsOnTrack())) {
348 rot = &crot->rioOnTrack(0);
355 const Surface *surf = &rot->associatedSurface();
356 const Amg::Vector3D measdir = surf->transform().rotation().col(0);
358 const double dotprod2 = measdir.dot(
360 surf->center().perp());
361 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
369 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
371 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
375 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
377 if (!cache.m_field_cache.solenoidOn()) {
378 const AmgVector(5) & newpars = parforcalo->parameters();
380 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
381 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
385 std::vector < MaterialEffectsOnTrack > calomeots;
388 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
392 parforcalo->associatedSurface(),
405 if (calomeots.empty()) {
410 std::unique_ptr<Track>
track;
413 cache.m_calomat =
false;
414 const bool tmp2 = cache.m_extmat;
415 const bool tmp4 = cache.m_idmat;
420 const double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
421 const double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
423 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
424 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
428 (errmatid !=
nullptr) &&
429 (errmatmuon !=
nullptr) &&
434 const double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
435 const double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
436 const double energyerror = std::sqrt(
437 calomeots[1].energyLoss()->sigmaDeltaE() *
438 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
439 pmuonerror * pmuonerror
443 (std::abs(calomeots[1].energyLoss()->deltaE()) -
444 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
447 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
452 parforcalo->associatedSurface(),
457 if (calomeots.empty()) {
464 const int nfits = cache.m_fit_status[
S_FITS];
465 bool firstfitwasattempted =
false;
468 if (!caloEntranceIsValid) {
473 (!cache.m_field_cache.toroidOn() && !cache.m_field_cache.solenoidOn()) ||
475 cache.m_getmaterialfromtrack &&
479 qoverpid * qoverpmuon > 0
484 if (cache.m_fit_status[
S_FITS] == (
unsigned int) (nfits + 1)) {
485 firstfitwasattempted =
true;
490 (
track ==
nullptr) &&
491 !firstfitwasattempted &&
492 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn())
495 GXFTrajectory trajectory2;
496 trajectory2.m_straightline = trajectory.m_straightline;
497 trajectory2.m_fieldprop = trajectory.m_fieldprop;
498 trajectory = trajectory2;
502 bool pseudoupdated =
false;
504 if (
track !=
nullptr) {
505 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
506 if (pseudostate ==
nullptr) {
517 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
522 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
524 pseudostate->measurement()->localParameters(),
525 pseudostate->measurement()->localCovariance()
528 if (updpar ==
nullptr) {
535 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
540 pseudopar->associatedSurface()
543 pseudostate->setMeasurement(std::move(newpseudo));
547 pseudostate->setMeasurementErrors(
errors);
548 pseudoupdated =
true;
552 trajectory.setConverged(
false);
553 cache.m_matfilled =
true;
559 *
track->perigeeParameters(),
561 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
564 cache.m_matfilled =
false;
568 cache.m_fit_status[
S_FITS] = nfits + 1;
570 if (
track !=
nullptr) {
571 track->info().addPatternReco(intrk1.info());
572 track->info().addPatternReco(intrk2.info());
576 cache.m_calomat =
tmp;
577 cache.m_extmat =
tmp2;
578 cache.m_idmat = tmp4;
◆ fitIm()
Definition at line 1849 of file GlobalChi2Fitter.cxx.
1857 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1859 GXFTrajectory trajectory;
1862 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1867 if (inputTrack.trackStateOnSurfaces()->empty()) {
1872 if (inputTrack.trackParameters()->empty()) {
1873 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1877 std::unique_ptr<const TrackParameters> minpar =
unique_clone(inputTrack.perigeeParameters());
1881 if (minpar ==
nullptr) {
1882 minpar =
unique_clone(*(inputTrack.trackParameters()->begin()));
1885 const bool tmpgetmat = cache.m_getmaterialfromtrack;
1891 cache.m_getmaterialfromtrack =
false;
1897 trajectory.trackStates().reserve(inputTrack.trackStateOnSurfaces()->size());
1899 const Surface *firsthitsurf =
nullptr;
1900 const Surface *lasthitsurf =
nullptr;
1902 bool hasmuon =
false;
1903 bool iscombined =
false;
1904 bool seenphimeas =
false;
1908 for (; itStates != endState; ++itStates) {
1909 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1911 (pMeasurement ==
nullptr) &&
1912 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1913 matEffects != inputTrack.info().particleHypothesis()
1918 if (pMeasurement !=
nullptr) {
1919 const Surface *surf = &pMeasurement->associatedSurface();
1920 Identifier hitid = surf->associatedDetectorElementIdentifier();
1924 hitid = crot->rioOnTrack(0).identify();
1928 if (firsthitsurf ==
nullptr) {
1929 firsthitsurf = surf;
1937 if ((**itStates).trackParameters() !=
nullptr) {
1938 lastidpar = (**itStates).trackParameters();
1939 if (firstidpar ==
nullptr) {
1940 firstidpar = lastidpar;
1945 const Amg::Vector3D measdir = surf->transform().rotation().col(0);
1946 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
1947 const double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1948 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1950 if (std::abs(surf->center().z()) > 13000) {
1953 if (surf->center().perp() > 9000 && std::abs(surf->center().z()) < 13000) {
1965 if (iscombined && seenphimeas && (phiem || phibo)) {
1975 cache.m_getmaterialfromtrack && trajectory.numberOfScatterers() != 0 &&
1976 (hasmuon || cache.m_acceleration)
1978 cache.m_matfilled =
true;
1981 if (firstidpar == lastidpar) {
1982 firstidpar = lastidpar =
nullptr;
1987 !cache.m_matfilled &&
1989 m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier()) !=
1992 (firstidpar !=
nullptr)
1994 if (
m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier())) {
2001 const bool tmpacc = cache.m_acceleration;
2003 const bool tmpsirecal = cache.m_sirecal;
2004 std::unique_ptr<Track> tmptrack =
nullptr;
2008 cache.m_fiteloss =
true;
2009 cache.m_sirecal =
false;
2012 cache.m_asymeloss =
true;
2015 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2016 cache.m_sirecal = tmpsirecal;
2018 if (tmptrack ==
nullptr) {
2019 cache.m_matfilled =
false;
2020 cache.m_getmaterialfromtrack = tmpgetmat;
2021 cache.m_acceleration = tmpacc;
2022 cache.m_fiteloss = tmpfiteloss;
2027 bool isbrem =
false;
2029 unsigned int n_brem=0;
2031 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
2032 GXFMaterialEffects *meff = state->materialEffects();
2034 if (meff !=
nullptr) {
2038 const MaterialProperties *matprop = meff->materialProperties();
2040 const double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->delta_p());
2042 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->position()));
2043 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2044 double costracksurf = 1.;
2046 costracksurf = std::abs(layerNormal.dot(layerpars->momentum().unit()));
2048 const double oldde = meff->deltaE();
2050 std::unique_ptr<EnergyLoss> eloss;
2051 double sigmascat = 0;
2053 if (matprop !=
nullptr) {
2054 eloss = std::make_unique<EnergyLoss>(
2055 m_elosstool->energyLoss(*matprop,
p, 1. / costracksurf,
2057 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop,
p, 1. / costracksurf, matEffects));
2059 if (eloss !=
nullptr) {
2060 meff->setDeltaE(eloss->deltaE());
2063 MaterialProperties
const tmpprop(1., meff->x0(), 0., 0., 0., 0.);
2064 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop,
p, 1. / costracksurf, matEffects));
2067 meff->setScatteringSigmas(
2075 meff->setDeltaE(oldde);
2076 if (!meff->isKink()) {
2077 meff->setSigmaDeltaE(0);
2080 bremdp = meff->delta_p();
2082 }
else if (eloss !=
nullptr) {
2083 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
2085 if (meff->sigmaDeltaE() > 0) {
2091 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2092 minpar=trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2093 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2097 cache.m_matfilled =
true;
2101 trajectory.brems().clear();
2103 trajectory.brems().resize(1);
2104 trajectory.brems()[0] = bremdp;
2107 cache.m_asymeloss =
false;
2108 trajectory.setNumberOfScatterers(nscats);
2111 trajectory.setNumberOfBrems(n_brem);
2115 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2117 bool pseudoupdated =
false;
2119 if ((
track !=
nullptr) && hasid && hasmuon) {
2120 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
2122 (pseudostate ==
nullptr) ||
2124 pseudostate->fitQuality().chiSquared() < 10
2130 const std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2132 pseudostate->measurement()->localParameters(),
2133 pseudostate->measurement()->localCovariance()
2136 if (updpar ==
nullptr) {
2143 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2146 pseudopar->associatedSurface()
2149 pseudostate->setMeasurement(std::move(newpseudo));
2153 pseudostate->setMeasurementErrors(
errors);
2154 pseudoupdated =
true;
2157 if (pseudoupdated) {
2158 trajectory.setConverged(
false);
2159 cache.m_matfilled =
true;
2161 cache.m_matfilled =
false;
2165 cache.m_matfilled =
false;
2166 cache.m_getmaterialfromtrack = tmpgetmat;
2167 cache.m_acceleration = tmpacc;
2168 cache.m_fiteloss = tmpfiteloss;
2170 if (
track !=
nullptr) {
2172 const TrackInfo& old_info = inputTrack.info();
2173 track->info().addPatternReco(old_info);
2176 return track.release();
◆ holesearchExtrapolation()
Helper method which performs an extrapolation with additional logic for hole search.
This method is a wrapper around extrapolateStepwise from the extrapolator interface, with the added functionality that it will null any returned track parameters which are on the start and end surface.
- Parameters
-
| [in] | ctx | An event context for extrapolation. |
| [in] | src | The track parameters to start extrapolating from. |
| [in] | dst | The track state to extrapolate to. |
| [in] | propdir | The propagation direction. |
- Returns
- A vector of track states, just like normal extrapolation.
Definition at line 7662 of file GlobalChi2Fitter.cxx.
7672 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7673 ctx,
src, dst.associatedSurface(), propdir,
false
7685 &
rv.front()->associatedSurface() == &dst.associatedSurface() ||
7686 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7687 trackParametersClose(*
rv.front(),
src, 0.001) ||
7688 trackParametersClose(*
rv.front(), *dst.trackParameters(), 0.001)
7691 rv.front().reset(
nullptr);
7701 &
rv.back()->associatedSurface() == &dst.associatedSurface() ||
7702 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7703 trackParametersClose(*
rv.back(),
src, 0.001) ||
7704 trackParametersClose(*
rv.back(), *dst.trackParameters(), 0.001)
7707 rv.back().reset(
nullptr);
◆ holeSearchHelper()
Helper method for the hole search that does the actual counting of holes and dead modules.
This is a helper function that does a lot of esoteric and weird things that you most likely won't need to know about. The gist of it is that you pass it a vector of track parameters and a counting object, and it will update those counters according to its analysis of the track parameters.
Unfortunately, due to the design of this method, it requires quite a lot of persistent state between invocations for the same track. That's bad design of course, but it is how it is for now. This means that there are quite a few state parameters.
- Parameters
-
| [in] | hc | A list of candidate hole track parameters to analyse. |
| [in,out] | id_set | A set of identifiers found to be holes or dead. |
| [in,out] | sct_set | A set of identifiers of SCT holes. |
| [in,out] | rv | The hole count container to update. |
| [in] | count_holes | Holes are counted only if this is enabled. |
| [in] | count_dead | Dead modules are counted only if this is enabled. |
Definition at line 7205 of file GlobalChi2Fitter.cxx.
7218 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7224 if (
tp ==
nullptr) {
7233 const TrkDetElementBase * de =
tp->associatedSurface().associatedDetectorElement();
7235 if (de ==
nullptr) {
7246 if (id_set.find(
id) != id_set.end()) {
7295 if (sct_set.find(
os) != sct_set.end()) {
7296 ++
rv.m_sct_double_hole;
◆ holeSearchProcess()
Conduct a hole search between a list of states, possibly reusing existing information.
Given a collection of state references, this method will conduct a hole search between consecutive pairs of states, possibly reusing existing information stored in the state data types. The method will check whether the state contains any previous hole search data and use it. If there is no data, it will run additional extrapolations to gather that data. It will then use a helper method to count holes and dead modules and return a total count.
In some cases, this method may error. Should this occur, it will return a non-extant value.
- Parameters
-
| [in] | ctx | An event context used for extrapolation. |
| [in] | states | A list of states to operate on, using consecutive states as extrapolation regions. |
- Returns
- A list of hole counts if the process succeeded, or a non-extant value in case of an error.
Definition at line 7390 of file GlobalChi2Fitter.cxx.
7404 constexpr
uint min_meas = 3;
7405 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7409 bool seen_meas =
false;
7411 std::set<Identifier> id_set;
7412 std::set<Identifier> sct_set;
7418 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7441 const double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7443 const bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7445 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7446 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7456 if (seen_meas && dist >= 2.5 && zStartValid) {
7463 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7464 std::vector<std::unique_ptr<TrackParameters>>
states;
7472 if (hc.has_value()) {
7493 GXFTrackState
const& last =
states.back();
7504 last.trackParameters() !=
nullptr &&
7511 std::vector<std::unique_ptr<Trk::TrackParameters>>
const bl =
m_extrapolator->extrapolateBlindly(
7513 *last.trackParameters(),
◆ holeSearchStates()
| std::vector< std::reference_wrapper< GXFTrackState > > Trk::GlobalChi2Fitter::holeSearchStates |
( |
GXFTrajectory & |
trajectory | ) |
const |
|
private |
Extracts a collection of track states which are important for hole search.
This method helps extract the measurement (and outlier) states from a track. These are the states between which we want to do a hole search, so the result of calling this method can be used as a source of truth for conducting a hole search on the track.
This method only returns states between the first and last measurements on the track, which is the region in which we are interested in doing a hole search.
As an example, if we denote scatterers as S, and measurements as M, this method would reduce the following track with numbered states:
1 2 3 4 5 6 7 8 9 M S S M M S M S S
Into a list of references [1, 4, 5, 7].
This method ensures that each pair of consecutive states in the return value list is a target for a hole search extrapolation.
- Parameters
-
| [in] | trajectory | The trajectory from which to extract states. |
- Returns
- A vector of state references as described above.
Definition at line 7315 of file GlobalChi2Fitter.cxx.
7323 GXFTrackState * lastmeas =
nullptr;
7325 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7337 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7338 rv.reserve(trajectory.trackStates().size());
7345 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7359 rv.emplace_back(*
s);
7366 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7368 if (de !=
nullptr) {
7381 if (
s.get() == lastmeas) {
◆ initFieldCache()
| void Trk::GlobalChi2Fitter::initFieldCache |
( |
const EventContext & |
ctx, |
|
|
Cache & |
cache |
|
) |
| const |
|
private |
Initialize a field cache inside a fit cache object.
Following the shift from old-style magnetic field services to the new cached implementation for thread safety, we need some additional logic to create a magnetic field cache object and insert it into our fitting cache object for access.
- Parameters
-
| [in] | cache | The GX2F cache objects in which to load the magnetic field cache. |
Definition at line 8510 of file GlobalChi2Fitter.cxx.
8520 if (cond_obj ==
nullptr) {
8525 cond_obj->getInitializedCache(cache.m_field_cache);
◆ initialize()
| StatusCode Trk::GlobalChi2Fitter::initialize |
( |
| ) |
|
|
overridevirtual |
Definition at line 210 of file GlobalChi2Fitter.cxx.
230 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
231 return StatusCode::FAILURE;
256 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
280 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
281 return StatusCode::FAILURE;
286 return StatusCode::SUCCESS;
◆ isMuonTrack()
| bool Trk::GlobalChi2Fitter::isMuonTrack |
( |
const Track & |
intrk1 | ) |
const |
|
private |
Definition at line 8468 of file GlobalChi2Fitter.cxx.
8469 const auto *pDataVector = intrk1.measurementsOnTrack();
8470 auto nmeas1 = pDataVector->size();
8471 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8479 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8481 if (lastMeasIsCompetingRIO){
8483 testrot = &testcrot->rioOnTrack(0);
8487 if (testrot ==
nullptr) {
8488 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8491 if(penultimateIsRIO){
8492 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8494 if (penultimateIsCompetingRIO){
8496 testrot = &testcrot->rioOnTrack(0);
8503 (testrot !=
nullptr) &&
◆ iterationsOfLastFit()
| int Trk::GlobalChi2Fitter::iterationsOfLastFit |
( |
| ) |
const |
|
privatevirtual |
◆ mainCombinationStrategy()
Definition at line 582 of file GlobalChi2Fitter.cxx.
590 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
595 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
596 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
598 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
599 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
600 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
602 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
606 if (muontrack->trackStateOnSurfaces()->empty()) {
612 muontrack->trackStateOnSurfaces()->end() - 1 :
613 muontrack->trackStateOnSurfaces()->begin();
615 const MeasurementBase *closestmuonmeas =
nullptr;
616 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
618 while (closestmuonmeas ==
nullptr) {
619 closestmuonmeas =
nullptr;
622 if ((**tsosit).measurementOnTrack() !=
nullptr) {
623 closestmuonmeas = (**tsosit).measurementOnTrack();
625 if (thispar !=
nullptr) {
626 const AmgVector(5) & parvec = thispar->parameters();
627 tp_closestmuon=thispar->associatedSurface().createUniqueTrackParameters(
628 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
642 std::unique_ptr<const TrackParameters> tmppar;
645 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
647 ctx, *tp_closestmuon, *cache.m_msEntrance, propdir,
nonInteracting);
650 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
652 if (tmppar !=
nullptr) {
653 const Surface & associatedSurface = tmppar->associatedSurface();
654 std::unique_ptr<Surface> muonsurf =
nullptr;
658 const CylinderBounds *cylbounds =
static_cast <const CylinderBounds *
>(&associatedSurface.bounds());
660 const double radius = cylbounds->r();
661 const double hlength = cylbounds->halflengthZ();
662 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
666 const double newz = (
667 associatedSurface.center().z() > 0 ?
668 associatedSurface.center().z() + 1 :
669 associatedSurface.center().z() - 1
673 associatedSurface.center().x(),
674 associatedSurface.center().y(),
678 trans.translation() << newpos;
680 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&associatedSurface.bounds());
681 const double rmin = discbounds->rMin();
682 const double rmax = discbounds->rMax();
683 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
687 if (muonsurf !=
nullptr) {
699 std::vector<const TrackStateOnSurface *> tmp_matvec;
701 if ((matvec !=
nullptr) && !matvec->empty()) {
702 tmp_matvec = *matvec;
703 delete tmp_matvec.back();
704 tmp_matvec.pop_back();
706 for (
auto &
i : tmp_matvec) {
711 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(
i->materialEffectsOnTrack());
713 const Surface *matsurf = &meff->associatedSurface();
720 trajectory.m_fieldprop,
725 if (tmppar ==
nullptr) {
733 trajectory.m_fieldprop,
739 if (tmppar ==
nullptr) {
747 const double de = std::abs(meff->energyLoss()->deltaE());
748 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
749 const double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
756 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
757 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
773 for (; itStates != endState; ++itStates) {
778 const bool tmpgetmat = cache.m_getmaterialfromtrack;
780 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
782 cache.m_extmat =
false;
784 cache.m_idmat =
false;
787 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
791 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
792 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
793 cache.m_getmaterialfromtrack =
true;
799 cache.m_getmaterialfromtrack = tmpgetmat;
806 trajectory.trackStates().back()->setTrackParameters(
nullptr);
809 std::unique_ptr<const TrackParameters> firstscatpar;
810 std::unique_ptr<const TrackParameters> lastscatpar;
813 double newqoverpid = 0;
816 const double de = std::abs(calomeots[1].energyLoss()->deltaE());
817 const double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
819 const double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
820 const double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
821 const double elosspull = (pbefore - pafter - de) / sigmade;
823 if (std::abs(elosspull) > 10) {
824 if (elosspull > 10) {
825 newqoverpid = 1 / (de + pafter + 10 * sigmade);
827 newqoverpid = 1 / (de + pafter - 10 * sigmade);
830 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
834 const AmgVector(5) & newpar = firstidpar->parameters();
835 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
836 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
844 if (lastidpar ==
nullptr) {
850 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
851 calomeots[0].associatedSurface(),
854 trajectory.m_fieldprop,
858 if (firstscatpar ==
nullptr) {
864 *(firstismuon ? firstidpar : tp_closestmuon),
865 calomeots[2].associatedSurface(),
868 trajectory.m_fieldprop,
872 if (lastscatpar ==
nullptr) {
876 std::optional<TransportJacobian> jac1;
877 std::optional<TransportJacobian> jac2;
878 std::unique_ptr<const TrackParameters> elosspar;
880 double firstscatphi = 0;
881 double secondscatphi = 0;
882 double firstscattheta = 0;
883 double secondscattheta = 0;
884 double muonscatphi = 0;
885 double muonscattheta = 0;
887 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
888 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
890 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
892 const Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
895 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
896 std::unique_ptr<const TrackParameters> startPar =
unique_clone(cache.m_idmat ? lastidpar.get() : indettrack->perigeeParameters());
898 for (
int i = 0;
i < 2;
i++) {
899 std::unique_ptr<const TrackParameters> tmpelosspar;
901 params1[
Trk::
phi] += muonscatphi;
902 params1[
Trk::
theta] += muonscattheta;
908 const std::unique_ptr<const TrackParameters> tmppar1(muonscatpar->associatedSurface().createUniqueTrackParameters(
909 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
921 trajectory.m_fieldprop,
930 calomeots[1].associatedSurface(),
932 trajectory.m_fieldprop
936 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
945 const AmgVector(5) & newpars = tmpelosspar->parameters();
946 const std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
947 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
951 elosspar = std::move(tmpelosspar);
954 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
958 calomeots[0].associatedSurface() :
959 calomeots[2].associatedSurface(),
962 trajectory.m_fieldprop,
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
977 trajectory.m_fieldprop
981 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
986 for (
int j = 0; j < 5; j++) {
987 for (
int k = 0;
k < 5;
k++) {
989 for (
int l = 0;
l < 5;
l++) {
990 jac3[j][
k] += (*jac2) (j,
l) * (*jac1) (
l,
k);
999 jac4(0, 0) = jac3[0][2];
1000 jac4(1, 1) = jac3[1][3];
1001 jac4(0, 1) = jac3[0][3];
1002 jac4(1, 0) = jac3[1][2];
1004 jac4 = jac4.inverse();
1016 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1018 if (cylsurf !=
nullptr) {
1022 if (discsurf !=
nullptr) {
1026 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1027 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1033 muonscatphi += dphi;
1034 muonscattheta += dtheta;
1037 const double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1040 firstscatphi = muonscatphi;
1041 secondscatphi = idscatphi;
1042 firstscattheta = muonscattheta;
1043 secondscattheta = idscattheta;
1045 firstscatphi = -idscatphi;
1046 secondscatphi = -muonscatphi;
1047 firstscattheta = -idscattheta;
1048 secondscattheta = -muonscattheta;
1051 if (
i == 1 && cache.m_field_cache.toroidOn() && !firstismuon) {
1053 params2[
Trk::
phi] += idscatphi;
1054 params2[
Trk::
theta] += idscattheta;
1060 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1061 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1063 idscatpar = firstscatpar.get();
1067 *cache.m_caloEntrance,
1071 if (startPar !=
nullptr) {
1080 rot.col(2) = trackdir;
1083 trans.linear().matrix() << rot;
1084 trans.translation() << startPar->position() - .1 * trackdir;
1085 PlaneSurface
const curvlinsurf(trans);
1095 if (curvlinpar !=
nullptr) {
1096 startPar.reset(curvlinpar);
1100 firstscatpar = std::move(scat2);
1104 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1105 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1106 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1108 const double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 const double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1112 for (
auto &
i : tmp_matvec) {
1117 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1118 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1121 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1123 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1126 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1128 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1129 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1130 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1133 for (
auto &
i : tmp_matvec) {
1140 if (startPar ==
nullptr) {
1145 (pull1 > 5 || pull2 > 5) &&
1151 bool largegap =
false;
1152 double previousz = 0;
1154 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1155 const MaterialEffectsBase *meff = (*itStates2)->materialEffectsOnTrack();
1157 const MeasurementBase *meas = (*itStates2)->measurementOnTrack();
1159 if (meff !=
nullptr) {
1161 const MaterialEffectsOnTrack *mefot{};
1163 mefot =
static_cast<const MaterialEffectsOnTrack *
>(meff);
1165 if ( mefot and mefot->energyLoss() and
1166 std::abs(mefot->energyLoss()->deltaE()) > 250 &&
1167 mefot->energyLoss()->sigmaDeltaE() < 1.e-9
1172 cache.m_extmat =
false;
1174 cache.m_idmat =
false;
1180 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1189 itStates2 == endState2 - 1 &&
1191 tpar->position().perp() > 9000 &&
1192 std::abs(tpar->position().z()) < 13000
1194 std::unique_ptr<const TrackParameters> pseudopar(tpar->clone());
1198 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1201 pseudopar->associatedSurface()
1204 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1211 pseudostate->setMeasurementErrors(
errors);
1212 trajectory.addMeasurementState(std::move(pseudostate));
1218 std::abs(trajectory.trackStates().back()->position().z()) > 20000 &&
1219 std::abs(previousz) < 12000
1225 previousz = trajectory.trackStates().back()->position().z();
1235 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
◆ makePerigee()
Definition at line 4471 of file GlobalChi2Fitter.cxx.
4476 const PerigeeSurface *persurf =
nullptr;
4479 persurf =
static_cast<const PerigeeSurface *
>(¶m.associatedSurface());
4481 if ((persurf !=
nullptr) && (!cache.m_acceleration || persurf->center().perp() > 5)) {
4483 return param.associatedSurface().createUniqueTrackParameters(
4488 if (cache.m_acceleration) {
4492 PerigeeSurface
const tmppersf;
4493 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4494 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4496 if (per ==
nullptr) {
4497 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4501 if(std::abs(per->position().z())>5000.) {
4502 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
◆ makeProtoState()
Definition at line 2557 of file GlobalChi2Fitter.cxx.
2569 ) && cache.m_getmaterialfromtrack
2571 if (cache.m_acceleration && trajectory.numberOfHits() == 0) {
2577 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(tsos->materialEffectsOnTrack());
2579 std::unique_ptr<GXFMaterialEffects> newmeff;
2582 meff->scatteringAngles() or
2583 meff->energyLoss() or
2585 (tsos->trackParameters() ==
nullptr)
2587 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2591 const double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2593 std::abs(1. / tsos->trackParameters()->parameters()[
Trk::qOverP]),
2606 meff->thicknessInX0(),
2612 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2616 (meff->energyLoss() !=
nullptr) &&
2617 meff->energyLoss()->sigmaDeltaE() > 0 &&
2620 ((meff->scatteringAngles() ==
nullptr) || meff->thicknessInX0() == 0)
2623 newmeff->setSigmaDeltaE(meff->energyLoss()->sigmaDeltaE());
2626 (tsos->trackParameters() !=
nullptr) &&
2627 !trajectory.trackStates().empty() &&
2628 ((**trajectory.trackStates().rbegin()).trackParameters() !=
nullptr)
2630 const double delta_p = 1000 * (
2631 tsos->trackParameters()->parameters()[
Trk::qOverP] -
2632 (**trajectory.trackStates().rbegin()).trackParameters()->
2636 newmeff->setdelta_p(delta_p);
2640 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(newmeff),
unique_clone(tsos->trackParameters())),
index);
2652 tsos->measurementOnTrack(),
2653 tsos->trackParameters(),
◆ makeProtoStateFromMeasurement()
Definition at line 2660 of file GlobalChi2Fitter.cxx.
2670 if (!measbase->associatedSurface().associatedDetectorElementIdentifier().is_valid()) {
2672 seg =
static_cast<const Segment *
>(measbase);
2679 imax = (
int) seg->numberOfMeasurementBases();
2682 for (
int i = 0;
i <
imax;
i++) {
2683 const MeasurementBase *measbase2 = ((seg !=
nullptr) &&
m_decomposesegments) ? seg->measurement(
i) : measbase;
2685 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2686 std::unique_ptr<const MeasurementBase>(measbase2->clone()),
2687 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->clone() :
nullptr)
2690 double sinstereo = 0;
2694 Identifier hitid = measbase2->associatedSurface().associatedDetectorElementIdentifier();
2699 hitid = crot->rioOnTrack(0).identify();
2703 bool measphi =
false;
2706 bool rotated =
false;
2722 if (measbase->localParameters().parameterKey() != 1) {
2743 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
2744 errors[0] = std::sqrt(covEigenValueSmall);
2745 sinstereo =
std::sin(covStereoAngle);
2758 const Surface *surf = &measbase2->associatedSurface();
2759 const Amg::Vector3D measdir = surf->transform().rotation().col(0);
2760 const double dotprod1 = measdir.dot(
Amg::Vector3D(0, 0, 1));
2761 const double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
2762 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2770 int param_index = 0;
2777 errors[1] = std::sqrt(
covmat(param_index, param_index));
2782 errors[2] = std::sqrt(
covmat(param_index, param_index));
2787 errors[3] = std::sqrt(
covmat(param_index, param_index));
2792 errors[4] = std::sqrt(
covmat(param_index, param_index));
2797 ATH_MSG_DEBUG(
"PseudoMeasurement, pos=" << measbase2->globalPosition());
2800 ATH_MSG_DEBUG(
"VertexOnTrack, pos=" << measbase2->globalPosition());
2803 ATH_MSG_DEBUG(
"Segment, pos=" << measbase2->globalPosition());
2813 ptsos->setMeasurementErrors(
errors);
2814 ptsos->setSinStereo(sinstereo);
2815 ptsos->setMeasurementType(hittype);
2816 ptsos->setMeasuresPhi(measphi);
2818 if (isoutlier && !cache.m_reintoutl) {
2823 const bool ok = trajectory.addMeasurementState(std::move(ptsos),
index);
◆ makeTrack()
Definition at line 7531 of file GlobalChi2Fitter.cxx.
7538 auto trajectory = std::make_unique<Trk::TrackStates>();
7544 GXFTrajectory tmptrajectory(oldtrajectory);
7546 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7548 if (perigee_ts ==
nullptr) {
7552 tmptrajectory.addBasicState(std::move(perigee_ts), cache.m_acceleration ? 0 : tmptrajectory.numberOfUpstreamStates());
7554 trajectory->reserve(tmptrajectory.trackStates().size());
7555 for (
auto & hit : tmptrajectory.trackStates()) {
7560 hit->resetTrackCovariance();
7566 hit->materialEffects()))
7572 auto trackState = hit->trackStateOnSurface();
7573 hit->resetTrackCovariance();
7574 trajectory->emplace_back(trackState.release());
7577 auto qual = std::make_unique<FitQuality>(tmptrajectory.chi2(), tmptrajectory.nDOF());
7588 if (matEffects ==
electron && tmptrajectory.hasKink()) {
7593 if (tmptrajectory.m_straightline) {
7597 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7607 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7616 std::optional<TrackHoleCount> hole_count;
7641 if (hole_count.has_value()) {
7654 rv->setTrackSummary(std::move(
ts));
◆ makeTrackFillDerivativeMatrix()
| void Trk::GlobalChi2Fitter::makeTrackFillDerivativeMatrix |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
oldtrajectory |
|
) |
| |
|
staticprivate |
Definition at line 6990 of file GlobalChi2Fitter.cxx.
6994 Amg::MatrixX & derivs = oldtrajectory.weightedResidualDerivatives();
6998 for (
auto & hit : oldtrajectory.trackStates()) {
6999 if (
const auto *pMeas{hit->measurement()};
7005 nrealmeas += hit->numberOfMeasuredParameters();
7008 cache.m_derivmat.resize(nrealmeas, oldtrajectory.numberOfFitParameters());
7009 cache.m_derivmat.setZero();
7012 const int nperpars = oldtrajectory.numberOfPerigeeParameters();
7013 const int nscat = oldtrajectory.numberOfScatterers();
7014 for (
auto & hit : oldtrajectory.trackStates()) {
7015 if (
const auto *pMeas{hit->measurement()};
7021 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
7022 for (
int j = 0; j < oldtrajectory.numberOfFitParameters(); j++) {
7023 cache.m_derivmat(
i, j) = derivs(measindex2, j) *
errors[measindex2];
7024 if ((j == 4 && !oldtrajectory.m_straightline) || j >= nperpars + 2 * nscat) {
7025 cache.m_derivmat(
i, j) *= 1000;
7032 measindex += hit->numberOfMeasuredParameters();
7033 }
else if (hit->materialEffects() ==
nullptr) {
7034 measindex2 += hit->numberOfMeasuredParameters();
◆ makeTrackFindPerigee()
◆ makeTrackFindPerigeeParameters()
Definition at line 7039 of file GlobalChi2Fitter.cxx.
7045 GXFTrackState *firstmeasstate =
nullptr;
7046 GXFTrackState *lastmeasstate =
nullptr;
7047 std::tie(firstmeasstate, lastmeasstate) = oldtrajectory.findFirstLastMeasurement();
7048 std::unique_ptr<const TrackParameters> per(
nullptr);
7051 std::unique_ptr<const TrackParameters> prevpar(
7052 firstmeasstate->trackParameters() !=
nullptr ?
7053 firstmeasstate->trackParameters()->clone() :
7056 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.upstreamMaterialLayers();
7060 if (prevpar ==
nullptr) {
7065 const Layer *
layer = layer1 !=
nullptr ? layer1 : layer2;
7067 const DistanceSolution distsol =
layer->surfaceRepresentation().straightLineDistanceEstimate(
7068 prevpar->position(), prevpar->momentum().unit()
7072 if (distsol.numberOfSolutions() == 2) {
7077 if (distsol.first() * distsol.second() < 0 && !
first) {
7086 std::unique_ptr<const TrackParameters> layerpar(
7090 layer->surfaceRepresentation(),
7093 oldtrajectory.m_fieldprop,
7098 if (layerpar ==
nullptr) {
7102 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7106 prevpar = std::move(layerpar);
7110 const Layer *startlayer = firstmeasstate->trackParameters()->associatedSurface().associatedLayer();
7112 if ((startlayer !=
nullptr) && (startlayer->layerMaterialProperties() !=
nullptr)) {
7113 double startfactor = startlayer->layerMaterialProperties()->alongPostFactor();
7114 const Surface & discsurf = startlayer->surfaceRepresentation();
7117 startfactor = startlayer->layerMaterialProperties()->oppositePostFactor();
7119 if (startfactor > 0.5) {
7120 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7121 firstmeasstate->trackParameters(), *startlayer,
oppositeMomentum, matEffects
7124 if (updatedpar !=
nullptr) {
7125 firstmeasstate->setTrackParameters(std::move(updatedpar));
7134 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
7136 if ((endlayer !=
nullptr) && (endlayer->layerMaterialProperties() !=
nullptr)) {
7137 double endfactor = endlayer->layerMaterialProperties()->alongPreFactor();
7138 const Surface & discsurf = endlayer->surfaceRepresentation();
7141 endfactor = endlayer->layerMaterialProperties()->oppositePreFactor();
7144 if (endfactor > 0.5) {
7145 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7146 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
7149 if (updatedpar !=
nullptr) {
7150 lastmeasstate->setTrackParameters(std::move(updatedpar));
7155 if (prevpar !=
nullptr) {
7162 oldtrajectory.m_fieldprop,
7167 if (per ==
nullptr) {
7168 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7173 }
else if (cache.m_acceleration && (firstmeasstate->trackParameters() !=
nullptr)) {
7175 *firstmeasstate->trackParameters(),
7181 per.reset(oldtrajectory.referenceParameters()->clone());
◆ myfit()
Definition at line 4509 of file GlobalChi2Fitter.cxx.
4517 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4521 if (!trajectory.m_straightline) {
4522 if (trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
4523 trajectory.m_straightline = !cache.m_field_cache.solenoidOn();
4524 }
else if ((trajectory.prefit() == 0) && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0) {
4525 trajectory.m_straightline = !cache.m_field_cache.toroidOn();
4527 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
4532 cache.m_lastiter = 0;
4536 if (trajectory.numberOfPerigeeParameters() == -1) {
4537 cache.incrementFitStatus(
S_FITS);
4538 if (trajectory.m_straightline) {
4539 trajectory.setNumberOfPerigeeParameters(4);
4541 trajectory.setNumberOfPerigeeParameters(5);
4545 if (trajectory.nDOF() < 0) {
4550 cache.m_phiweight.clear();
4551 cache.m_firstmeasurement.clear();
4552 cache.m_lastmeasurement.clear();
4555 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4559 if (matEffects ==
Trk::electron && trajectory.m_straightline) {
4565 trajectory.setMass(
mass);
4567 ATH_MSG_DEBUG(
"start param: " << param <<
" pos: " << param.position() <<
" pt: " << param.pT());
4569 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4571 if (!cache.m_acceleration && (per ==
nullptr)) {
4581 cache.m_acceleration &&
4582 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits() &&
4585 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4588 cache.m_fastmat =
false;
4593 !cache.m_acceleration ||
4594 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() != trajectory.numberOfHits()
4596 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4599 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4603 if (cache.m_acceleration && (trajectory.referenceParameters() ==
nullptr) && (per ==
nullptr)) {
4606 if (trajectory.numberOfScatterers() >= 2) {
4607 GXFTrackState *scatstate =
nullptr;
4608 GXFTrackState *scatstate2 =
nullptr;
4611 for (
const auto & state : trajectory.trackStates()) {
4614 scatindex == trajectory.numberOfScatterers() / 2 ||
4615 state->materialEffects()->deltaE() == 0
4617 scatstate2 = state.get();
4622 scatstate = state.get();
4629 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4630 throw std::logic_error(
"Invalid scatterer");
4633 vertex = .49 * (scatstate->position() + scatstate2->position());
4635 const int nstates = (
int) trajectory.trackStates().size();
4637 trajectory.trackStates()[nstates / 2 - 1]->position() +
4638 trajectory.trackStates()[nstates / 2]->position()
4642 PerigeeSurface
const persurf(
vertex);
4643 std::unique_ptr<const TrackParameters> nearestpar;
4644 double mindist = 99999;
4645 std::vector < GXFTrackState * >mymatvec;
4647 for (
auto &
it : trajectory.trackStates()) {
4648 if ((*it).trackParameters() ==
nullptr) {
4653 .straightLineDistanceEstimate(
4654 (*it).trackParameters()->position(),
4655 (*it).trackParameters()->momentum().unit())
4658 const bool insideid = (
4659 (cache.m_caloEntrance ==
nullptr) ||
4660 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4664 (((*it).measurement() !=
nullptr) && insideid) || (
4665 ((*it).materialEffects() !=
nullptr) &&
4667 (*it).materialEffects()->deltaE() == 0 ||
4668 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4670 (*it).materialEffects()->deltaPhi() != 0
4674 const double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4675 if (dist < mindist) {
4683 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4684 mymatvec.push_back(
it.get());
4688 if (nearestpar ==
nullptr) {
4692 for (
auto & state : mymatvec) {
4694 const Surface &matsurf = state->associatedSurface();
4695 const DistanceSolution distsol = matsurf.straightLineDistanceEstimate(
4696 nearestpar->position(), nearestpar->momentum().unit());
4700 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4704 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4710 trajectory.m_fieldprop,
4714 if (tmppar ==
nullptr) {
4722 trajectory.m_fieldprop,
4726 if (tmppar ==
nullptr) {
4738 if (state->materialEffects()->sigmaDeltaE() > 0) {
4739 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4742 const double de = std::abs(state->materialEffects()->deltaE());
4743 const double oldp = std::abs(1 / newpars[
Trk::qOverP]);
4744 const double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4750 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4751 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4755 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4761 trajectory.m_fieldprop,
4766 if (tmpPars !=
nullptr) {
4767 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4772 const double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4773 const double toteloss = std::abs(trajectory.totalEnergyLoss());
4774 const double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4778 per = per->associatedSurface().createUniqueTrackParameters(
4783 if (per ==
nullptr) {
4791 PerigeeSurface
const persurf2(per->position());
4792 per = persurf2.createUniqueTrackParameters(
4800 }
else if (per ==
nullptr) {
4804 if ((per ==
nullptr) && (trajectory.referenceParameters() ==
nullptr)) {
4812 if (trajectory.m_straightline && (per !=
nullptr)) {
4813 if (trajectory.numberOfPerigeeParameters() == -1) {
4814 trajectory.setNumberOfPerigeeParameters(4);
4818 per = per->associatedSurface().createUniqueTrackParameters(
4821 }
else if (trajectory.numberOfPerigeeParameters() == -1) {
4822 trajectory.setNumberOfPerigeeParameters(5);
4825 if (per !=
nullptr) {
4826 trajectory.setReferenceParameters(std::move(per));
4829 const int nfitpar = trajectory.numberOfFitParameters();
4830 const int nperpars = trajectory.numberOfPerigeeParameters();
4831 const int nscat = trajectory.numberOfScatterers();
4832 const int nbrem = trajectory.numberOfBrems();
4835 Eigen::MatrixXd a_inv;
4836 a.resize(nfitpar, nfitpar);
4841 derivPool.setZero();
4843 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
4844 if (state->materialEffects() !=
nullptr) {
4847 state->setDerivatives(derivPool);
4850 bool doderiv =
true;
4851 const int tmpminiter = cache.m_miniter;
4854 cache.m_lastiter =
it;
4860 cache.m_miniter = tmpminiter;
4864 if (!trajectory.converged()) {
4865 cache.m_fittercode =
4875 cache.m_miniter = tmpminiter;
4879 const int nhits = trajectory.numberOfHits();
4880 const int ntrthits = trajectory.numberOfTRTHits();
4881 const int nsihits = trajectory.numberOfSiliconHits();
4882 const double redchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
4883 const double prevredchi2 = (trajectory.nDOF() > 0) ? trajectory.prevchi2() / trajectory.nDOF() : 0;
4892 (redchi2 < prevredchi2 &&
4893 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4894 nsihits + ntrthits == nhits
4899 if (
it != 1 || nsihits != 0 || trajectory.nDOF() <= 0 || trajectory.chi2() / trajectory.nDOF() <= 3) {
4904 cache.m_miniter = tmpminiter;
4911 const int ntrtprechits = trajectory.numberOfTRTPrecHits();
4912 const int ntrttubehits = trajectory.numberOfTRTTubeHits();
4914 if (ntrtprechits+ntrttubehits) {
4915 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
4917 if (phf<m_minphfcut && it>=3) {
4918 if ((ntrtprechits+ntrttubehits)>=15) {
4923 <<
" | nTRTPrecHits = " << ntrtprechits
4924 <<
" | nTRTTubeHits = " << ntrttubehits
4925 <<
" | nOutliers = "
4926 << trajectory.numberOfOutliers());
4928 if (!trajectory.converged()) {
4934 cache.m_miniter = tmpminiter;
4943 cache.m_miniter = tmpminiter;
4945 if (trajectory.prefit() == 0) {
4949 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
4950 if (lltOfW.info() == Eigen::Success) {
4954 const int ncols =
a.cols();
4955 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4956 a_inv = lltOfW.solve(weightInvAMG);
4965 GXFTrajectory *finaltrajectory = &trajectory;
4967 (runOutlier || cache.m_sirecal) &&
4968 trajectory.numberOfSiliconHits() == trajectory.numberOfHits()
4975 if (traj != &trajectory) {
4980 finaltrajectory = traj;
4987 if (!cache.m_acceleration && (finaltrajectory->prefit() == 0)) {
4988 if (nperpars == 5) {
4989 for (
int i = 0;
i <
a.cols();
i++) {
4990 a_inv(4,
i) *= .001;
4991 a_inv(
i, 4) *= .001;
4995 int scatterPos = nperpars + 2 * nscat;
4996 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4997 for (
int i = 0;
i <
a.cols();
i++) {
4998 a_inv(scatterPos,
i) *= .001;
4999 a_inv(
i, scatterPos) *= .001;
5005 const int nperparams = finaltrajectory->numberOfPerigeeParameters();
5006 for (
int i = 0;
i < nperparams;
i++) {
5007 for (
int j = 0; j < nperparams; j++) {
5008 (errmat) (j,
i) = a_inv(j,
i);
5012 if (trajectory.m_straightline) {
5013 (errmat) (4, 4) = 1
e-20;
5016 const AmgVector(5) & perpars = finaltrajectory->referenceParameters()->parameters();
5017 std::unique_ptr<const TrackParameters> measper(
5018 finaltrajectory->referenceParameters()->associatedSurface().createUniqueTrackParameters(
5019 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5023 finaltrajectory->setReferenceParameters(std::move(measper));
5025 cache.m_fullcovmat = a_inv;
5029 std::unique_ptr<Track>
track =
nullptr;
5031 if (finaltrajectory->prefit() > 0) {
5032 if (finaltrajectory != &trajectory) {
5034 delete finaltrajectory;
5039 if (finaltrajectory->numberOfOutliers() <=
m_maxoutliers || !runOutlier) {
5046 const double cut = (finaltrajectory->numberOfSiliconHits() ==
5047 finaltrajectory->numberOfHits())
5053 (
track !=
nullptr) && (
5054 track->fitQuality()->numberDoF() != 0 &&
5055 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5058 track.reset(
nullptr);
5062 if (
track ==
nullptr) {
5066 if (finaltrajectory != &trajectory) {
5067 delete finaltrajectory;
5070 return track.release();
◆ myfit_helper()
◆ numericalDerivatives()
Definition at line 8305 of file GlobalChi2Fitter.cxx.
8319 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8322 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8329 const Surface & previousSurface = tmpprevpar->associatedSurface();
8333 for (
int i = 0;
i < 5;
i++) {
8336 if (thisdiscsurf &&
i == 1) {
8342 if (
i == 0 && thiscylsurf) {
8344 }
else if (
i == 1 && thisdiscsurf) {
8350 std::unique_ptr<const TrackParameters> parpluseps(
8351 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8360 const std::unique_ptr<const TrackParameters> parminuseps(
8361 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8371 std::unique_ptr<const TrackParameters> newparpluseps(
8382 std::unique_ptr<const TrackParameters> newparminuseps(
8397 if (newparpluseps ==
nullptr) {
8409 if (newparminuseps ==
nullptr) {
8421 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8425 for (
int j = 0; j < 5; j++) {
8429 if (j == 0 && cylsurf) {
8431 }
else if (j == 1 && discsurf) {
8435 (*jac) (j,
i) =
diff / (2 * eps[
i]);
◆ processTrkVolume()
Definition at line 2833 of file GlobalChi2Fitter.cxx.
2837 if (tvol ==
nullptr) {
2844 if (confinedLayers !=
nullptr) {
2848 if (
layer !=
nullptr) {
2853 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2857 const CylinderLayer *cyllay =
nullptr;
2859 cyllay =
static_cast<const CylinderLayer *
>(
layer);
2862 const DiscLayer *disclay =
nullptr;
2864 disclay =
static_cast<const DiscLayer *
>(
layer);
2867 if (disclay !=
nullptr) {
2868 if (disclay->center().z() < 0) {
2869 cache.m_negdiscs.push_back(disclay);
2871 cache.m_posdiscs.push_back(disclay);
2873 }
else if (cyllay !=
nullptr) {
2874 cache.m_barrelcylinders.push_back(cyllay);
2884 for (
size_t ib = 0 ;
ib < bsurf.size(); ++
ib) {
2885 const Layer *
layer = bsurf[
ib]->surfaceRepresentation().materialLayer();
2887 if (
layer ==
nullptr)
continue;
2891 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2895 const CylinderSurface *cylsurf =
nullptr;
2898 cylsurf =
static_cast<const CylinderSurface *
>(&
layer->surfaceRepresentation());
2900 const DiscSurface *discsurf =
nullptr;
2903 discsurf =
static_cast<const DiscSurface *
>(&
layer->surfaceRepresentation());
2905 if (discsurf !=
nullptr) {
2907 discsurf->center().z() < 0 &&
2908 std::find(cache.m_negdiscs.begin(), cache.m_negdiscs.end(),
layer) == cache.m_negdiscs.end()
2910 cache.m_negdiscs.push_back(
layer);
2912 discsurf->center().z() > 0 &&
2913 std::find(cache.m_posdiscs.begin(), cache.m_posdiscs.end(),
layer) == cache.m_posdiscs.end()
2915 cache.m_posdiscs.push_back(
layer);
2918 (cylsurf !=
nullptr) &&
2919 std::find(cache.m_barrelcylinders.begin(), cache.m_barrelcylinders.end(),
layer) == cache.m_barrelcylinders.end()
2921 cache.m_barrelcylinders.push_back(
layer);
2924 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2931 if (confinedVolumes !=
nullptr) {
2932 for (
const auto & volume : confinedVolumes->arrayObjects()) {
2933 if (volume !=
nullptr) {
◆ retrieveTrackingGeometry()
◆ runIteration()
Definition at line 5994 of file GlobalChi2Fitter.cxx.
6004 const int nDOFold = trajectory.nDOF();
6005 const double oldChi2 = trajectory.chi2();
6006 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6008 if (cache.m_phiweight.empty()) {
6009 cache.m_phiweight.assign(trajectory.trackStates().size(), 1);
6028 int bremno_maxbrempull = 0;
6029 GXFTrackState* state_maxbrempull =
nullptr;
6044 if ((state_maxbrempull !=
nullptr) && trajectory.converged()) {
6045 trajectory.setConverged(
false);
6046 trajectory.setChi2(1e15);
6053 const int nDOFnew = trajectory.nDOF();
6054 const double newChi2 = trajectory.chi2();
6055 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6057 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6058 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6060 if (trajectory.prefit() > 0 && trajectory.converged()) {
6069 if (cache.m_firstmeasurement.empty()) {
6073 if (
a.cols() != trajectory.numberOfFitParameters()) {
6093 if (doDeriv || weightChanged) {
6104 if (trajectory.prefit() == 0) {
6105 if (trajectory.converged()) {
6106 const int nSiHits = trajectory.numberOfSiliconHits();
6107 const int nTrtHits = trajectory.numberOfTRTHits();
6108 const int nHits = trajectory.numberOfHits();
6110 if (
nSiHits + nTrtHits != nHits) {
6117 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
◆ runTrackCleanerSilicon()
- Warning
- This method has some unclear memory ownership mechanics that might not correspond fully with the model described at the beginning of the file. Be aware!
Definition at line 6511 of file GlobalChi2Fitter.cxx.
6520 bool trackok =
false;
6521 GXFTrajectory *oldtrajectory = &trajectory;
6522 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6523 GXFTrajectory *newtrajectory =
nullptr;
6524 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6530 while (!trackok && oldtrajectory->nDOF() > 0) {
6532 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->trackStates();
6535 Amg::MatrixX & weightderiv = oldtrajectory->weightedResidualDerivatives();
6536 const int nfitpars = oldtrajectory->numberOfFitParameters();
6537 const int nhits = oldtrajectory->numberOfHits();
6538 const int nsihits = oldtrajectory->numberOfSiliconHits();
6540 if (nhits != nsihits) {
6544 double maxsipull = -1;
6546 int hitno_maxsipull = -1;
6547 int measno_maxsipull = -1;
6548 int stateno_maxsipull = 0;
6549 GXFTrackState *state_maxsipull =
nullptr;
6554 const int noutl = oldtrajectory->numberOfOutliers();
6560 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6561 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6567 double *
errors = state->measurementErrors();
6569 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6570 const double sinstereo = state->sinStereo();
6571 const double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6572 double weight1 = -1;
6574 if (hitcov(0, 0) > trackcov(0, 0)) {
6575 if (sinstereo == 0) {
6579 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6580 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6585 const double weight2 = (
6591 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6592 const double sipull2 = (
6594 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6597 sipull1 =
std::max(sipull1, sipull2);
6599 if (sipull1 > maxsipull) {
6600 maxsipull = sipull1;
6601 measno_maxsipull = measno;
6602 state_maxsipull = state.get();
6603 stateno_maxsipull = stateno;
6604 hitno_maxsipull = hitno;
6615 measno += state->numberOfMeasuredParameters();
6619 const double maxpull = maxsipull;
6621 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6622 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6631 oldtrajectory->chi2() / oldtrajectory->nDOF() > .25 *
m_chi2cut
6633 state_maxsipull = oldtrajectory->trackStates()[stateno_maxsipull].get();
6634 const PrepRawData *prd{};
6636 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6637 prd = rot->prepRawData();
6639 std::unique_ptr < const RIO_OnTrack > broadrot;
6640 double *olderror = state_maxsipull->measurementErrors();
6642 const TrackParameters *trackpar_maxsipull = state_maxsipull->trackParameters();
6644 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6645 const std::unique_ptr<const TrackParameters> trackparForCorrect(
6646 trackpar_maxsipull->associatedSurface().createUniqueTrackParameters(
6652 state_maxsipull->hasTrackCovariance()
6654 state_maxsipull->trackCovariance())
6658 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6659 double newpull = -1;
6660 double newpull1 = -1;
6661 double newpull2 = -1;
6662 double newres1 = -1;
6663 double newres2 = -1;
6664 double newsinstereo = 0;
6668 !state_maxsipull->isRecalibrated() &&
6670 oldtrajectory->chi2() / trajectory.nDOF() > .3 *
m_chi2cut &&
6679 if (state_maxsipull->sinStereo() != 0) {
6680 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
6681 newerror[0] = std::sqrt(covEigenValueSmall);
6682 newsinstereo =
std::sin(covStereoAngle);
6684 newerror[0] = std::sqrt(
covmat(0, 0));
6687 const double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6689 if (cosstereo != 1.) {
6691 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6692 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6695 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6698 if (newerror[0] == 0.0) {
6699 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6702 newpull1 = std::abs(newres1 / newerror[0]);
6706 newerror[1] = std::sqrt(
covmat(1, 1));
6707 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6708 newpull2 = std::abs(newres2 / newerror[1]);
6711 newpull =
std::max(newpull1, newpull2);
6717 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6719 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6720 throw std::runtime_error(
6721 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6726 newtrajectory = oldtrajectory;
6728 if (
a.cols() != nfitpars) {
6732 const double oldres1 =
res[measno_maxsipull];
6733 res[measno_maxsipull] = newres1;
6734 err[measno_maxsipull] = newerror[0];
6736 for (
int i = 0;
i < nfitpars;
i++) {
6737 if (weightderiv(measno_maxsipull,
i) == 0) {
6741 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6743 for (
int j =
i; j < nfitpars; j++) {
6747 weightderiv(measno_maxsipull,
i) *
6748 weightderiv(measno_maxsipull, j) *
6749 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6753 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6757 const double oldres2 =
res[measno_maxsipull + 1];
6758 res[measno_maxsipull + 1] = newres2;
6759 err[measno_maxsipull + 1] = newerror[1];
6761 for (
int i = 0;
i < nfitpars;
i++) {
6762 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6766 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6768 for (
int j =
i; j < nfitpars; j++) {
6772 weightderiv(measno_maxsipull + 1,
i) *
6773 weightderiv(measno_maxsipull + 1, j) *
6774 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6779 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6784 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6785 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6786 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6787 olderror[1] <<
" newerror_1=" << newerror[1]
6790 state_maxsipull->setMeasurement(std::move(broadrot));
6791 state_maxsipull->setSinStereo(newsinstereo);
6792 state_maxsipull->setMeasurementErrors(newerror);
6796 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6797 (oldtrajectory->chi2() / oldtrajectory->nDOF() > .3 *
m_chi2cut || noutl > 1)
6801 (oldtrajectory->nDOF() > 1 || hittype_maxsipull ==
TrackState::SCT) &&
6806 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6807 measno_maxsipull <<
" pull=" << maxsipull
6814 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6815 newtrajectory = cleanup_newtrajectory.get();
6817 if (newa.cols() != nfitpars) {
6822 Amg::MatrixX & newweightderiv = newtrajectory->weightedResidualDerivatives();
6823 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6824 throw std::runtime_error(
6825 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6829 const double oldres1 =
res[measno_maxsipull];
6830 newres[measno_maxsipull] = 0;
6832 for (
int i = 0;
i < nfitpars;
i++) {
6833 if (weightderiv(measno_maxsipull,
i) == 0) {
6837 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6839 for (
int j =
i; j < nfitpars; j++) {
6843 weightderiv(measno_maxsipull,
i) *
6844 weightderiv(measno_maxsipull, j)
6848 newweightderiv(measno_maxsipull,
i) = 0;
6852 const double oldres2 =
res[measno_maxsipull + 1];
6853 newres[measno_maxsipull + 1] = 0;
6855 for (
int i = 0;
i < nfitpars;
i++) {
6856 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6860 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6862 for (
int j =
i; j < nfitpars; j++) {
6863 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6870 weightderiv(measno_maxsipull + 1,
i) *
6871 weightderiv(measno_maxsipull + 1, j)
6875 newweightderiv(measno_maxsipull + 1,
i) = 0;
6879 newtrajectory->setOutlier(stateno_maxsipull);
6885 newtrajectory->setConverged(
false);
6901 if (!newtrajectory->converged()) {
6903 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6910 if (!newtrajectory->converged()) {
6919 const double oldchi2 = oldtrajectory->chi2() / oldtrajectory->nDOF();
6920 const double newchi2 = (newtrajectory->nDOF() > 0) ? newtrajectory->chi2() / newtrajectory->nDOF() : 0;
6923 if (newtrajectory->nDOF() != oldtrajectory->nDOF() && maxsipull > cut2) {
6924 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6926 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6931 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6932 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6940 (void)cleanup_oldtrajectory.release();
6941 return oldtrajectory;
6943 if (oldtrajectory != newtrajectory) {
6944 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6945 oldtrajectory = newtrajectory;
6952 Eigen::LLT < Eigen::MatrixXd >
const lltOfW(
a);
6953 if (lltOfW.info() == Eigen::Success) {
6957 const int ncols =
a.cols();
6958 Amg::MatrixX const weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6959 fullcov = lltOfW.solve(weightInvAMG);
6977 oldtrajectory->nDOF() > 0 &&
6978 oldtrajectory->chi2() / oldtrajectory->nDOF() >
m_chi2cut &&
6986 (void)cleanup_oldtrajectory.release();
6987 return oldtrajectory;
◆ runTrackCleanerTRT()
Definition at line 6341 of file GlobalChi2Fitter.cxx.
6354 if (
it == 1 && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
6358 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6361 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6362 const int nfitpars = trajectory.numberOfFitParameters();
6364 if (
a.cols() != nfitpars) {
6368 const int nperpars = trajectory.numberOfPerigeeParameters();
6369 const int nscats = trajectory.numberOfScatterers();
6372 bool outlierremoved =
false;
6373 bool hitrecalibrated =
false;
6375 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6376 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6384 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6388 trajectory.setOutlier(stateno);
6389 outlierremoved =
true;
6391 double *
errors = state->measurementErrors();
6392 const double olderror =
errors[0];
6394 trajectory.updateTRTHitCount(stateno, olderror);
6396 for (
int i = 0;
i < nfitpars;
i++) {
6397 if (weightderiv(measno,
i) == 0) {
6401 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6403 for (
int j =
i; j < nfitpars; j++) {
6406 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6409 weightderiv(measno,
i) = 0;
6413 }
else if (trtrecal) {
6414 double *
errors = state->measurementErrors();
6415 const double olderror =
errors[0];
6417 const auto *
const thisMeasurement{state->measurement()};
6423 if (oldrot->prepRawData() !=
nullptr) {
6424 const double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6426 const double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6428 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6429 const double distance = std::abs(std::abs(trackradius) - dcradius);
6431 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6432 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6433 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6434 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6437 if (newrot !=
nullptr) {
6439 hitrecalibrated =
true;
6443 if ((measno < 0) or (measno >= (
int)
res.size())) {
6444 throw std::runtime_error(
6445 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6449 const double oldres =
res[measno];
6450 const double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6452 state->setMeasurement(std::move(newrot));
6454 trajectory.updateTRTHitCount(stateno, olderror);
6456 for (
int i = 0;
i < nfitpars;
i++) {
6457 if (weightderiv(measno,
i) == 0) {
6461 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6463 for (
int j =
i; j < nfitpars; j++) {
6467 !cache.m_phiweight.empty() &&
6470 i < nperpars + 2 * nscats &&
6471 (
i - nperpars) % 2 == 0
6473 weight = cache.m_phiweight[(
i - nperpars) / 2];
6478 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6481 weightderiv(measno,
i) *= olderror / newerror;
6484 res[measno] = newres;
6485 err[measno] = newerror;
6494 measno += state->numberOfMeasuredParameters();
6498 if (trajectory.nDOF() < 0) {
6503 if (outlierremoved || hitrecalibrated) {
6505 trajectory.setConverged(
false);
6507 cache.m_miniter =
it + 2;
◆ setMinIterations()
| void Trk::GlobalChi2Fitter::setMinIterations |
( |
int |
| ) |
|
|
privatevirtual |
Definition at line 8446 of file GlobalChi2Fitter.cxx.
8448 (
"Configure the minimum number of Iterations via jobOptions");
◆ throwFailedToGetTrackingGeomtry()
| void Trk::GlobalChi2Fitter::throwFailedToGetTrackingGeomtry |
( |
| ) |
const |
|
private |
◆ trackingGeometry()
Definition at line 1132 of file GlobalChi2Fitter.h.
1135 if (!cache.m_trackingGeometry)
1137 return cache.m_trackingGeometry;
◆ tryToConverge()
Definition at line 5404 of file GlobalChi2Fitter.cxx.
5411 const double oldChi2 = trajectory.prevchi2();
5412 const double newChi2 = trajectory.chi2();
5417 const double nDOF = trajectory.nDOF();
5418 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5419 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5422 trajectory.prefit() > 0 && (
5423 (newRedChi2 < 2 &&
it != 0) ||
5424 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 &&
it != 1)
5427 trajectory.setConverged(
true);
5433 const int nsihits = trajectory.numberOfSiliconHits();
5434 const int ntrthits = trajectory.numberOfTRTHits();
5435 const int nhits = trajectory.numberOfHits();
5437 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5438 miniter =
std::max(miniter, cache.m_miniter);
5440 if (
it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5441 trajectory.setConverged(
true);
◆ tryToWeightAfromMaterial()
Definition at line 5821 of file GlobalChi2Fitter.cxx.
5830 const int nPerPars = trajectory.numberOfPerigeeParameters();
5836 bool weightChanged =
false;
5843 double newPhiWeight = 1.1;
5844 double newThetaWeight = 1.001;
5845 if (trajectory.prefit() == 0) {
5851 newPhiWeight = 1.00000001;
5852 }
else if (
it == 1) {
5853 newPhiWeight = 1.0000001;
5854 }
else if (
it <= 3) {
5855 newPhiWeight = 1.0001;
5856 }
else if (
it <= 6) {
5857 newPhiWeight = 1.01;
5860 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5861 newPhiWeight = 1.0001;
5862 newThetaWeight = 1.0001;
5863 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5864 newPhiWeight = 1.001;
5865 newThetaWeight = 1.0001;
5872 std::size_t scatno = 0;
5877 for (
const auto & state : trajectory.trackStates()) {
5878 const GXFMaterialEffects *meff = state->materialEffects();
5880 if (meff ==
nullptr) {
5884 const bool isValidPlaneSurface =
5886 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5891 if (meff->deltaE() == 0 || (trajectory.prefit() == 0 && isValidPlaneSurface)) {
5892 weightChanged =
true;
5894 const int scatNoIndex = 2 * scatno + nPerPars;
5896 if (trajectory.prefit() == 0 && meff->sigmaDeltaPhi() != 0) {
5897 if (scatno >= cache.m_phiweight.size()) {
5899 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5900 throw std::range_error(
message.str());
5908 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5911 cache.m_phiweight[scatno] = newPhiWeight;
5912 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5913 }
else if (trajectory.prefit() >= 2) {
5914 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5915 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5927 meff->sigmaDeltaPhi() != 0 &&
5928 (trajectory.prefit() == 0 || meff->deltaE() == 0)
5942 trajectory.prefit() == 2 &&
5944 trajectory.numberOfBrems() > 0 &&
5945 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5950 return weightChanged;
◆ updateEnergyLoss()
Definition at line 7960 of file GlobalChi2Fitter.cxx.
7973 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7977 double newqoverp = 0;
7979 if (meff.sigmaDeltaE() <= 0) {
7984 const double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.deltaE()) * std::sqrt(
mass *
mass + oldp * oldp) + meff.deltaE() * meff.deltaE();
7991 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7997 return surf.createUniqueTrackParameters(
7998 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
◆ updateFitParameters()
Method to update peregee parameters, scattering angles, and brems.
Tries to solve the system [A] * deltaParameters = b and then update in the trajectory all parameters used for the fit. Returns also a status.
Definition at line 6126 of file GlobalChi2Fitter.cxx.
6139 Eigen::LLT<Eigen::MatrixXd>
const llt(lu_m);
6141 if (llt.info() != Eigen::Success) {
6151 const int nscat = trajectory.numberOfScatterers();
6152 const int nbrem = trajectory.numberOfBrems();
6153 const int nperparams = trajectory.numberOfPerigeeParameters();
6165 double d0 = refpar->parameters()[
Trk::d0];
6166 double z0 = refpar->parameters()[
Trk::z0];
6169 double qoverp = refpar->parameters()[
Trk::qOverP];
6171 if (nperparams > 0) {
6172 d0 += deltaParameters[0];
6173 z0 += deltaParameters[1];
6174 phi += deltaParameters[2];
6175 theta += deltaParameters[3];
6176 qoverp = (trajectory.m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6188 std::vector < std::pair < double, double >>&scatangles = trajectory.scatteringAngles();
6189 for (
int i = 0;
i < nscat;
i++) {
6190 scatangles[
i].first += deltaParameters[2 *
i + nperparams];
6191 scatangles[
i].second += deltaParameters[2 *
i + nperparams + 1];
6197 std::vector < double >&delta_ps = trajectory.brems();
6198 for (
int i = 0;
i < nbrem;
i++) {
6199 delta_ps[
i] += deltaParameters[nperparams + 2 * nscat +
i];
6205 std::unique_ptr<const TrackParameters> newper(
6206 trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
6214 trajectory.setReferenceParameters(std::move(newper));
6215 trajectory.setScatteringAngles(scatangles);
6216 trajectory.setBrems(delta_ps);
◆ updatePixelROTs()
Update the Pixel ROT using the current trajectory/local track parameters.
Definition at line 6221 of file GlobalChi2Fitter.cxx.
6227 if ( trajectory.numberOfSiliconHits() == 0) {
6236 if (!splitProbContainer.isValid()) {
6240 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6243 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6244 const int nfitpars = trajectory.numberOfFitParameters();
6247 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6252 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6255 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6265 const PrepRawData *prd{};
6267 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6268 prd = rot->prepRawData();
6278 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6279 if (!splitProb.isSplit()) {
6280 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6284 std::unique_ptr < const RIO_OnTrack > newrot;
6285 double *olderror = state->measurementErrors();
6288 double newerror[5] = {-1,-1,-1,-1,-1};
6289 double newres[2] = {-1,-1};
6291 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6298 newerror[0] = std::sqrt(
covmat(0, 0));
6299 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6300 newerror[1] = std::sqrt(
covmat(1, 1));
6301 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6303 if (
a.cols() != nfitpars) {
6308 for(
int k =0;
k<2;
k++ ){
6309 const double oldres =
res[measno+
k];
6310 res[measno+
k] = newres[
k];
6311 err[measno+
k] = newerror[
k];
6313 for (
int i = 0;
i < nfitpars;
i++) {
6314 if (weightderiv(measno+
k,
i) == 0) {
6318 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6320 for (
int j =
i; j < nfitpars; j++) {
6324 weightderiv(measno+
k,
i) *
6325 weightderiv(measno+
k, j) *
6326 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6330 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6334 state->setMeasurement(std::move(newrot));
6335 state->setMeasurementErrors(newerror);
◆ updateSystemWithMaxBremPull()
Definition at line 5445 of file GlobalChi2Fitter.cxx.
5453 if (state_maxbrempull ==
nullptr) {
5457 state_maxbrempull->materialEffects()->setSigmaDeltaE(
5458 10 * state_maxbrempull->materialEffects()->sigmaDeltaEPos()
5461 state_maxbrempull->materialEffects()->setKink(
true);
5463 const int nbrem = trajectory.numberOfBrems();
5465 const int nmeas = (
int)
res.size();
5468 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5469 const double newError = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5470 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5472 const int nFitPars = trajectory.numberOfFitParameters();
5473 if (
a.cols() != nFitPars) {
5477 const double errorRatio = oldError / newError;
5478 const double errorReductionRatio = 1 -
std::pow(errorRatio, 2);
5480 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5481 for (
int i = 0;
i < nFitPars;
i++) {
5482 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5486 for (
int j =
i; j < nFitPars; j++) {
5487 const double newaij =
a(
i, j) - errorReductionRatio *
5488 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5491 a.fillSymmetric(
i, j, newaij);
5493 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= errorRatio;
◆ ATLAS_THREAD_SAFE
| std::array<std::atomic<unsigned int>, S_MAX_VALUE> m_fit_status Trk::GlobalChi2Fitter::ATLAS_THREAD_SAFE = {} |
|
mutableprivate |
◆ m_acceleration
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_acceleration {this, "Acceleration", false} |
|
private |
◆ m_asymeloss
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_asymeloss {this, "AsymmetricEnergyLoss", true} |
|
private |
◆ m_boundaryCheckTool
| ToolHandle<IBoundaryCheckTool> Trk::GlobalChi2Fitter::m_boundaryCheckTool {this, "BoundaryCheckTool", "", "Boundary checking tool for detector sensitivities" } |
|
private |
◆ m_broadROTcreator
| ToolHandle<IRIO_OnTrackCreator> Trk::GlobalChi2Fitter::m_broadROTcreator {this, "BroadRotCreatorTool", "", ""} |
|
private |
◆ m_calomat
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_calomat {this, "MuidMat", false} |
|
private |
◆ m_caloMaterialProvider
◆ m_calotool
◆ m_calotoolparam
◆ m_chi2cut
| Gaudi::Property<double> Trk::GlobalChi2Fitter::m_chi2cut {this, "TrackChi2PerNDFCut", 1.e15} |
|
private |
◆ m_clusterSplitProbContainer
◆ m_createSummary
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_createSummary {this, "CreateTrackSummary", true} |
|
private |
◆ m_decomposesegments
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_decomposesegments {this, "DecomposeSegments", true} |
|
private |
◆ m_DetID
◆ m_domeastrackpar
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_domeastrackpar {this, "MeasuredTrackParameters", true} |
|
private |
◆ m_elosstool
◆ m_extensioncuts
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_extensioncuts {this, "TRTExtensionCuts", true} |
|
private |
◆ m_extmat
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_extmat {this, "ExtrapolatorMaterial", true} |
|
private |
◆ m_extrapolator
◆ m_field_cache_key
Initial value:{
this,
"AtlasFieldCacheCondObj",
"fieldCondObj",
"Trk::GlobalChi2Fitter field conditions object key"
}
Definition at line 1157 of file GlobalChi2Fitter.h.
◆ m_fillderivmatrix
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_fillderivmatrix {this, "FillDerivativeMatrix", false} |
|
private |
◆ m_fiteloss
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_fiteloss {this, "FitEnergyLoss", false} |
|
private |
◆ m_fixbrem
| Gaudi::Property<int> Trk::GlobalChi2Fitter::m_fixbrem {this, "FixBrem", -1} |
|
private |
◆ m_getmaterialfromtrack
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_getmaterialfromtrack {this, "GetMaterialFromTrack", true} |
|
private |
◆ m_holeSearch
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_holeSearch {this, "DoHoleSearch", false} |
|
private |
◆ m_idVolume
◆ m_kinkfinding
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_kinkfinding {this, "KinkFinding", false} |
|
private |
◆ m_matupdator
◆ m_maxit
| Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxit {this, "MaxIterations", 30} |
|
private |
◆ m_maxitPixelROT
| Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxitPixelROT {this, "IterationsToRebuildPixelRots", 0} |
|
private |
◆ m_maxoutliers
| Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxoutliers {this, "MaxOutliers", 10} |
|
private |
◆ m_miniter
| Gaudi::Property<int> Trk::GlobalChi2Fitter::m_miniter {this, "MinimumIterations", 1} |
|
private |
◆ m_minphfcut
| Gaudi::Property<double> Trk::GlobalChi2Fitter::m_minphfcut {this, "MinPHFCut", 0.} |
|
private |
◆ m_navigator
| ToolHandle<INavigator> Trk::GlobalChi2Fitter::m_navigator {this, "NavigatorTool", "Trk::Navigator/CosmicsNavigator", ""} |
|
private |
◆ m_numderiv
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_numderiv {this, "NumericalDerivs", false} |
|
private |
◆ m_outlcut
| Gaudi::Property<double> Trk::GlobalChi2Fitter::m_outlcut {this, "OutlierCut", 5.0} |
|
private |
◆ m_p
| Gaudi::Property<double> Trk::GlobalChi2Fitter::m_p {this, "Momentum", 0.0} |
|
private |
◆ m_propagator
| ToolHandle<IPropagator> Trk::GlobalChi2Fitter::m_propagator {this, "PropagatorTool", "", ""} |
|
private |
◆ m_redoderivs
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_redoderivs {this, "RecalculateDerivatives", false} |
|
private |
◆ m_reintoutl
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_reintoutl {this, "ReintegrateOutliers", false} |
|
private |
◆ m_rejectLargeNScat
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_rejectLargeNScat {this, "RejectLargeNScat", false} |
|
private |
◆ m_residualPullCalculator
◆ m_ROTcreator
| ToolHandle<IRIO_OnTrackCreator> Trk::GlobalChi2Fitter::m_ROTcreator {this, "RotCreatorTool", "", ""} |
|
private |
◆ m_scalefactor
| Gaudi::Property<double> Trk::GlobalChi2Fitter::m_scalefactor {this, "TRTTubeHitCut", 2.5} |
|
private |
◆ m_scattool
◆ m_signedradius
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_signedradius {this, "SignedDriftRadius", true} |
|
private |
◆ m_sirecal
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_sirecal {this, "RecalibrateSilicon", false} |
|
private |
◆ m_storemat
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_storemat {this, "StoreMaterialOnTrack", true} |
|
private |
◆ m_straightlineprop
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_straightlineprop {this, "StraightLine", true} |
|
private |
◆ m_trackingGeometryReadKey
Initial value:{
this,
"TrackingGeometryReadKey",
"AtlasTrackingGeometry",
"Key of the TrackingGeometry conditions data."
}
Definition at line 1150 of file GlobalChi2Fitter.h.
◆ m_trtrecal
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_trtrecal {this, "RecalibrateTRT", false} |
|
private |
◆ m_updator
| ToolHandle<IUpdator> Trk::GlobalChi2Fitter::m_updator {this, "MeasurementUpdateTool", "", ""} |
|
private |
◆ m_useCaloTG
| Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_useCaloTG {this, "UseCaloTG", false} |
|
private |
The documentation for this class was generated from the following files:
def retrieve(aClass, aKey=None)
bool is_pixel(Identifier id) const
std::unique_ptr< const TrackParameters > makePerigee(Cache &, const TrackParameters &, const ParticleHypothesis) const
void makeProtoState(Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const TrackingGeometry * retrieveTrackingGeometry(const EventContext &ctx) const
bool is_rpc(Identifier id) const
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Gaudi::Property< bool > m_createSummary
Const iterator class for DataVector/DataList.
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
ToolHandle< INavigator > m_navigator
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
void fillResidualsAndErrors(const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const
@ DeadElement
outside the element
Gaudi::Property< double > m_outlcut
@ z
global position (cartesian)
bool ensureValidEntranceMuonSpectrometer(const EventContext &ctx, Cache &cache) const
Gaudi::Property< bool > m_numderiv
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotool
std::string find(const std::string &s)
return a remapped string
ToolHandle< IMaterialEffectsUpdator > m_matupdator
bool is_csc(Identifier id) const
Scalar perp() const
perp method - perpenticular length
bool inside(const Amg::Vector3D &gp, double tol=0.) const
Inside() method for checks.
Gaudi::Property< bool > m_holeSearch
bool is_sct(Identifier id) const
const Amg::Vector3D & position() const
Access method for the position.
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
FitterStatusCode runIteration(const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const
@ numberOfSCTDeadSensors
number of TRT hits
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
bool isMuonTrack(const Track &) const
@ BremFitSuccessful
A brem fit was performed on this track and this fit was successful.
Gaudi::Property< bool > m_redoderivs
std::pair< long int, long int > indices
Gaudi::Property< bool > m_straightlineprop
static void fillAfromScatterers(GXFTrajectory &trajectory, Amg::SymMatrixX &a)
bool contains(ParamDefs par) const
The simple check for the clients whether the parameter is contained.
bool is_mm(Identifier id) const
PropagationResult calculateTrackParametersPropagate(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possib...
static void fillAfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
ToolHandle< IRIO_OnTrackCreator > m_ROTcreator
static bool tryToWeightAfromMaterial(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2)
ToolHandle< IExtrapolator > m_extrapolator
Track * myfit(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const
Gaudi::Property< bool > m_decomposesegments
represents a deflection of the track caused through multiple scattering in material.
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
ToolHandle< IPropagator > m_propagator
@ loc2
generic first and second local coordinate
void runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
void addMaterialUpdateTrajectory(Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * >> &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const
Given layer information, probe those layers for scatterers and add them to a track.
std::vector< size_t > vec
bool is_trt(Identifier id) const
std::vector< std::unique_ptr< TrackParameters > > holesearchExtrapolation(const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const
Helper method which performs an extrapolation with additional logic for hole search.
#define ATH_MSG_VERBOSE(x)
bool const RAWDATA *ch2 const
ToolHandle< IRIO_OnTrackCreator > m_broadROTcreator
base class to integrate material effects on Trk::Track in a flexible way.
std::unique_ptr< const TrackParameters > makeTrackFindPerigeeParameters(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
bool empty() const
Test if the key is blank.
Gaudi::Property< int > m_maxoutliers
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ numberOfSCTHoles
number of Holes in both sides of a SCT module
bool is_valid() const
Check if id is in a valid state.
SG::ReadCondHandleKey< TrackingGeometry > m_trackingGeometryReadKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_field_cache_key
ToolHandle< IBoundaryCheckTool > m_boundaryCheckTool
@ u
Enums for curvilinear frames.
FitterStatusCode calculateTrackParameters(const EventContext &ctx, GXFTrajectory &, bool) const
void stable_sort(std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, std::reverse_iterator< DataModel_detail::iterator< DVL > > end, Compare comp)
Specialization of stable_sort for DataVector/List.
Gaudi::Property< int > m_fixbrem
static void fillBfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b)
std::pair< double, ParamDefs > DefinedParameter
void initFieldCache(const EventContext &ctx, Cache &cache) const
Initialize a field cache inside a fit cache object.
@ MATERIAL_EFFECTS_ON_TRACK
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
void updateSystemWithMaxBremPull(GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
const LayerArray * confinedLayers() const
Return the subLayer array.
Gaudi::Property< int > m_maxit
@ StraightTrack
A straight track.
AmgSymMatrix(5) &GXFTrackState
@ SlimmedTrack
A slimmed track.
represents the full description of deflection and e-loss of a track in material.
ToolHandle< IResidualPullCalculator > m_residualPullCalculator
std::vector< std::reference_wrapper< GXFTrackState > > holeSearchStates(GXFTrajectory &trajectory) const
Extracts a collection of track states which are important for hole search.
static void addMaterialGetLayers(Cache &cache, std::vector< std::pair< const Layer *, const Layer * >> &layers, std::vector< std::pair< const Layer *, const Layer * >> &uplayers, const std::vector< std::unique_ptr< GXFTrackState >> &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat)
Collect all possible layers that a given track could have passed through.
static void makeTrackFillDerivativeMatrix(Cache &, GXFTrajectory &)
ToolHandle< Trk::ITrkMaterialProviderTool > m_caloMaterialProvider
double getDistance(const xAOD::Vertex *vtx1, const xAOD::Vertex *vtx2)
ToolHandle< IEnergyLossUpdator > m_elosstool
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
Track * mainCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
bool is_tgc(Identifier id) const
Gaudi::Property< bool > m_domeastrackpar
@ numberOfPixelDeadSensors
number of pixel hits with broad errors (width/sqrt(12))
Eigen::Affine3d Transform3D
const AtlasDetectorID * m_DetID
PropagationResult calculateTrackParametersPropagateHelper(const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const
Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method...
std::pair< std::vector< unsigned int >, bool > res
double chi2(TH1 *h0, TH1 *h1)
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionDisc(Cache &cache, const DiscSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a disc surface.
cut
This script demonstrates how to call a C++ class from Python Also how to use PyROOT is shown.
MeasurementType
enum describing the flavour of MeasurementBase
@ ExtrapolationFailure
extrapolation failed
void throwFailedToGetTrackingGeomtry() const
Gaudi::Property< bool > m_fiteloss
std::optional< GlobalChi2Fitter::TrackHoleCount > holeSearchProcess(const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState >> &states) const
Conduct a hole search between a list of states, possibly reusing existing information.
static bool correctAngles(double &, double &)
static std::optional< std::pair< Amg::Vector3D, double > > addMaterialFindIntersectionCyl(Cache &cache, const CylinderSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat)
Find the intersection of a set of track parameters onto a cylindrical surface.
Gaudi::Property< bool > m_extensioncuts
void holeSearchHelper(const std::vector< std::unique_ptr< TrackParameters >> &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const
Helper method for the hole search that does the actual counting of holes and dead modules.
const TrackingGeometry * trackingGeometry(Cache &cache, const EventContext &ctx) const
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotoolparam
BinnedArray< TrackingVolume > TrackingVolumeArray
void calculateTrackErrors(GXFTrajectory &, Amg::SymMatrixX &, bool) const
std::unique_ptr< GXFTrackState > makeTrackFindPerigee(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
@ BremFit
A brem fit was performed on this track.
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
@ ExtrapolationFailureDueToSmallMomentum
extrapolation failed due to small momentum
int value() const
layerIndex expressed in an integer
Gaudi::Property< bool > m_fillderivmatrix
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
void tryToConverge(const Cache &cache, GXFTrajectory &trajectory, const int it) const
ToolHandle< IMultipleScatteringUpdator > m_scattool
std::unique_ptr< Track > makeTrack(const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const
@ FullField
Field is set to be realistic, but within a given Volume.
bool processTrkVolume(Cache &, const Trk::TrackingVolume *tvol) const
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
@ BremPoint
This represents a brem point on the track, and so will contain TrackParameters and MaterialEffectsBas...
std::optional< TransportJacobian > numericalDerivatives(const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const
@ numberOfPixelHoles
number of pixels which have a ganged ambiguity.
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
Gaudi::Property< bool > m_acceleration
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
@ NoField
Field is set to 0., 0., 0.,.
Ensure that the ATLAS eigen extensions are properly loaded.
void addMaterial(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const
void updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const
Update the Pixel ROT using the current trajectory/local track parameters.
std::vector< std::string > remainder(const std::vector< std::string > &v1, const std::vector< std::string > &v2)
std::string to_string(const DetectorType &type)
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
GXFTrajectory * runTrackCleanerSilicon(const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const
ParametersBase< TrackParametersDim, Charged > TrackParameters
bool is_indet(Identifier id) const
double charge(const T &p)
Gaudi::Property< double > m_chi2cut
StatusCode initialize(bool used=true)
Eigen::Matrix< double, 3, 1 > Vector3D
FitterStatusCode updateFitParameters(GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const
Method to update peregee parameters, scattering angles, and brems.
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
static void compensatePhiWeights(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a)
const Amg::Vector3D & momentum() const
Access method for the momentum.
Track * backupCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
@ Biased
RP with track state including the hit.
bool is_muon(Identifier id) const
@ GlobalChi2Fitter
Track's from Thijs' global chi^2 fitter.
void makeProtoStateFromMeasurement(Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const
Gaudi::Property< double > m_p
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
#define ATH_MSG_WARNING(x)
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
def delta_phi(phi1, phi2)
Gaudi::Property< int > m_maxitPixelROT
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
Track * fitIm(const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const
std::unique_ptr< T > unique_clone(const T *v)
Gaudi::Property< double > m_scalefactor
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
ToolHandle< IUpdator > m_updator
void addIDMaterialFast(const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final
@ PseudoMeasurementOnTrack
bool is_stgc(Identifier id) const
Gaudi::Property< bool > m_calomat
Gaudi::Property< bool > m_rejectLargeNScat
static void fillFirstLastMeasurement(Cache &cache, GXFTrajectory &trajectory)
Gaudi::Property< bool > m_trtrecal
static constexpr std::array< ParamDefs, 6 > pardef
Constructor.
bool consistentSurfaces(U)
virtual MaterialEffectsDerivedType derivedType() const =0
Returns the concrete derived type.
constexpr int pow(int base, int exp) noexcept
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
bool ensureValidEntranceCalo(const EventContext &ctx, Cache &cache) const
Scalar mag() const
mag method
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > updateEnergyLoss(const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
static void calculateDerivatives(GXFTrajectory &)
float nSiHits(const U &p)
bool is_mdt(Identifier id) const
Gaudi::Property< bool > m_useCaloTG
void fillDerivatives(GXFTrajectory &traj) const
virtual double r() const override final
This method returns the radius.
@ OutlierLogicFailure
outlier logic failed
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_clusterSplitProbContainer
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
@ Unknown
Track fitter not defined.