![Logo](../../ATLAS-Logo-Square-Blue-RGB.png) |
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 | fillResiduals (const EventContext &ctx, Cache &, GXFTrajectory &, int, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool &) const |
|
void | fillDerivatives (GXFTrajectory &traj, bool onlybrem=false) const |
|
FitterStatusCode | runIteration (const EventContext &ctx, Cache &, GXFTrajectory &, int, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool &) const |
|
FitterStatusCode | updateFitParameters (GXFTrajectory &, Amg::VectorX &, const Amg::SymMatrixX &) const |
|
void | updatePixelROTs (GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &) 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 |
|
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 |
|
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 | 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 3532 of file GlobalChi2Fitter.cxx.
3543 if (cache.m_caloEntrance ==
nullptr) {
3547 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
3552 if (cache.m_caloEntrance ==
nullptr) {
3563 cache.m_negdiscs.empty() &&
3564 cache.m_posdiscs.empty() &&
3565 cache.m_barrelcylinders.empty()
3580 cache.m_fastmat =
false;
3581 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3597 bool hasmat =
false;
3598 int indexoffset = 0, lastmatindex = 0;
3599 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.trackStates();
3601 GXFTrackState *firstsistate =
nullptr;
3602 GXFTrackState *lastsistate =
nullptr;
3613 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3614 if (oldstates[
i]->materialEffects() !=
nullptr) {
3623 if (firstsistate ==
nullptr) {
3624 if (oldstates[
i]->trackParameters() ==
nullptr) {
3625 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3628 oldstates[
i]->associatedSurface(),
3631 trajectory.m_fieldprop,
3635 if (tmppar ==
nullptr)
return;
3637 oldstates[
i]->setTrackParameters(std::move(tmppar));
3639 firstsistate = oldstates[
i].get();
3641 lastsistate = oldstates[
i].get();
3649 if (lastsistate ==
nullptr) {
3650 throw std::logic_error(
"No track state");
3659 if (lastsistate->trackParameters() ==
nullptr) {
3660 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3663 lastsistate->associatedSurface(),
3665 trajectory.m_fieldprop,
3669 if (tmppar ==
nullptr)
return;
3671 lastsistate->setTrackParameters(std::move(tmppar));
3680 refpar = lastsistate->trackParameters();
3681 indexoffset = lastmatindex;
3683 refpar = firstsistate->trackParameters();
3696 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3697 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.upstreamMaterialLayers();
◆ addMaterial()
Definition at line 3710 of file GlobalChi2Fitter.cxx.
3717 if (refpar2 ==
nullptr) {
3720 const MeasurementBase *firstmuonhit =
nullptr;
3721 const MeasurementBase *lastmuonhit =
nullptr;
3722 const MeasurementBase *firstidhit =
3724 const MeasurementBase *lastidhit =
nullptr;
3725 const MeasurementBase *firsthit =
nullptr;
3726 const MeasurementBase *lasthit =
nullptr;
3727 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3728 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3729 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3730 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3731 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3732 bool matvec_used=
false;
3733 std::unique_ptr<TrackParameters> startmatpar1;
3734 std::unique_ptr<TrackParameters> startmatpar2;
3745 int npseudomuon1 = 0;
3746 int npseudomuon2 = 0;
3748 for (
auto & state :
states) {
3751 GXFMaterialEffects *meff = state->materialEffects();
3754 if (firstidhit ==
nullptr) {
3763 if (firsthit ==
nullptr) {
3764 firsthit = state->measurement();
3765 if (cache.m_acceleration) {
3766 if (
tp ==
nullptr) {
3770 state->associatedSurface(),
3776 if (
tp ==
nullptr) {
3780 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3786 lasthit = state->measurement();
3792 if (firstidhit ==
nullptr) {
3793 firstidhit = state->measurement();
3796 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3800 lastidhit = state->measurement();
3801 if (
tp !=
nullptr) {
3806 if (firstsiliconpar ==
nullptr) {
3807 firstsiliconpar =
tp;
3809 lastsiliconpar =
tp;
3821 if (firstmuonhit ==
nullptr) {
3822 firstmuonhit = state->measurement();
3823 if (
tp !=
nullptr) {
3827 lastmuonhit = state->measurement();
3828 if (
tp !=
nullptr) {
3834 if (meff->deltaE() == 0) {
3835 if (firstcalopar ==
nullptr) {
3836 firstcalopar = state->trackParameters();
3838 lastcalopar = state->trackParameters();
3840 if (firstmatpar ==
nullptr) {
3841 firstmatpar = state->trackParameters();
3846 std::unique_ptr<TrackParameters> refpar;
3849 if (trajectory.m_straightline &&
m_p != 0) {
3853 refpar = refpar2->associatedSurface().createUniqueTrackParameters(
3854 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3857 if (firstmatpar !=
nullptr) {
3862 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3866 double mass = trajectory.mass();
3868 const AmgVector(5) & newpars = startmatpar2->parameters();
3872 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3873 newpars[0], newpars[1], newpars[2], newpars[3],
3874 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3878 }
else if (trajectory.m_straightline &&
m_p != 0) {
3882 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3883 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3889 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3890 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3894 if ((firstidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3896 DistanceSolution distsol = firstidhit->associatedSurface().straightLineDistanceEstimate(
3898 refpar->momentum().unit()
3903 if (
distance < 0 && distsol.numberOfSolutions() > 0 && !cache.m_acceleration) {
3904 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3906 const Surface *destsurf = &firstidhit->associatedSurface();
3907 std::unique_ptr<const TrackParameters> tmppar;
3909 if (firstmuonhit !=
nullptr) {
3910 if (cache.m_caloEntrance ==
nullptr) {
3914 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
3920 if (cache.m_caloEntrance ==
nullptr) {
3925 *cache.m_caloEntrance,
3929 if (tmppar !=
nullptr) {
3930 destsurf = &tmppar->associatedSurface();
3935 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
3940 false, matEffects) );
3943 if (matvec && !matvec->empty()) {
3944 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3945 const MaterialEffectsBase *meb = (*matvec)[
i]->materialEffectsOnTrack();
3948 const MaterialEffectsOnTrack *meot =
static_cast < const MaterialEffectsOnTrack *
>(meb);
3949 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3950 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->clone() :
nullptr;
3951 meff->setSigmaDeltaE(0);
3952 matstates.push_back(std::make_unique<GXFTrackState>(
3954 std::unique_ptr<const TrackParameters>(newpars)
3964 if ((lastidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3965 DistanceSolution distsol = lastidhit->associatedSurface().straightLineDistanceEstimate(
3967 refpar->momentum().unit()
3972 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
3973 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3974 const Surface *destsurf = &lastidhit->associatedSurface();
3975 std::unique_ptr<const TrackParameters> tmppar;
3976 std::unique_ptr<Surface> calosurf;
3977 if (firstmuonhit !=
nullptr) {
3978 if (cache.m_caloEntrance ==
nullptr) {
3982 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
3988 if (cache.m_caloEntrance ==
nullptr) {
3993 *cache.m_caloEntrance,
3998 if (tmppar !=
nullptr) {
3999 const CylinderSurface *cylcalosurf =
nullptr;
4002 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
4004 const DiscSurface *disccalosurf =
nullptr;
4007 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
4009 if (cylcalosurf !=
nullptr) {
4011 const CylinderBounds & cylbounds = cylcalosurf->bounds();
4012 double radius = cylbounds.r();
4013 double hlength = cylbounds.halflengthZ();
4014 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
4015 }
else if (disccalosurf !=
nullptr) {
4017 disccalosurf->center().z() > 0 ?
4018 disccalosurf->center().z() - 1 :
4019 disccalosurf->center().z() + 1
4023 disccalosurf->center().x(),
4024 disccalosurf->center().y(),
4029 trans.translation() << newpos;
4031 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&disccalosurf->bounds());
4032 double rmin = discbounds->rMin();
4033 double rmax = discbounds->rMax();
4034 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
4036 destsurf = calosurf.release();
4040 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4042 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
4043 matvec_used =
false;
4045 if (matvec && !matvec->empty()) {
4046 for (
const auto &
i : *matvec) {
4051 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4052 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4053 if (cache.m_fiteloss && (meot->energyLoss() !=
nullptr)) {
4054 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4057 if (matEffects ==
electron && cache.m_asymeloss) {
4058 meff->setDeltaE(-5);
4060 if (trajectory.numberOfTRTHits() == 0) {
4061 meff->setScatteringSigmas(0, 0);
4064 meff->setSigmaDeltaE(50);
4067 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
4069 matstates.push_back(std::make_unique<GXFTrackState>(
4071 std::unique_ptr<const TrackParameters>(newparams)
4083 if (cache.m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
4086 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
4090 firstmuonhit->associatedSurface(),
4095 if (calomeots.empty()) {
4100 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4101 if (lasthit == lastmuonhit) {
4102 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4105 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4108 calomeots[
i].associatedSurface(),
4111 trajectory.m_fieldprop,
4115 if (layerpar ==
nullptr) {
4120 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4123 lastcalopar = layerpar.get();
4127 double qoverp = layerpar->parameters()[
Trk::qOverP];
4128 double qoverpbrem = 0;
4132 (firstmuonpar !=
nullptr) &&
4133 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4135 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4137 double sign = (qoverp > 0) ? 1 : -1;
4138 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4141 const AmgVector(5) & newpar = layerpar->parameters();
4143 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4144 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4146 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4149 matstates.push_back(std::make_unique<GXFTrackState>(
4151 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4153 prevtrackpars = std::move(layerpar);
4158 firsthit == firstmuonhit &&
4159 (!cache.m_getmaterialfromtrack || lasthit == lastidhit)
4162 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4164 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4167 calomeots[
i].associatedSurface(),
4170 trajectory.m_fieldprop,
4174 if (layerpar ==
nullptr) {
4179 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4188 double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4193 (lastmuonpar !=
nullptr) &&
4194 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4198 double sign = (qoverpbrem > 0) ? 1 : -1;
4199 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4202 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4203 const AmgVector(5) & newpar = layerpar->parameters();
4205 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4206 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4210 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4215 if (lasthit == lastmuonhit && cache.m_extmat) {
4216 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4218 if (lastcalopar !=
nullptr) {
4219 if (cache.m_msEntrance ==
nullptr) {
4223 cache.m_msEntrance =
geometry->trackingVolume(
"MuonSpectrometerEntrance");
4229 if (cache.m_msEntrance ==
nullptr) {
4231 }
else if (cache.m_msEntrance->inside(lastcalopar->position())) {
4234 *cache.m_msEntrance,
4238 if (muonpar1 !=
nullptr) {
4246 rot.col(2) = trackdir;
4248 trans.linear().matrix() << rot;
4249 trans.translation() << muonpar1->
position() - .1 * trackdir;
4250 PlaneSurface curvlinsurf(trans);
4252 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4260 if (curvlinpar !=
nullptr) {
4261 muonpar1 = std::move(curvlinpar);
4265 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->clone());
4268 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4271 DistanceSolution distsol;
4273 if (muonpar1 !=
nullptr) {
4274 distsol = lastmuonhit->associatedSurface().straightLineDistanceEstimate(
4282 if ((
distance > 0) and(distsol.numberOfSolutions() >
4283 0) and (firstmuonhit !=
nullptr)) {
4284 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4291 if (distsol.numberOfSolutions() == 1) {
4293 }
else if (distsol.numberOfSolutions() == 2) {
4295 std::abs(distsol.first()) < std::abs(distsol.second()) ?
4301 if (
distance < 0 && distsol.numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4302 if (firstmuonpar !=
nullptr) {
4305 if (trajectory.m_straightline &&
m_p != 0) {
4310 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4313 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4316 firstmuonhit->associatedSurface(),
4319 trajectory.m_fieldprop,
4323 if (tmppar !=
nullptr) {
4324 muonpar1 = std::move(tmppar);
4330 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4331 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4334 states.back()->associatedSurface(),
4338 matvec_used =
false;
4345 if (matvec && !matvec->empty()) {
4346 for (
int j = 0; j < (
int) matvec->size(); j++) {
4347 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4351 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4352 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4355 !trajectory.m_straightline &&
4356 (meot->energyLoss() !=
nullptr) &&
4357 std::abs(meff->deltaE()) > 25 &&
4358 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4360 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4363 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4365 matstates.push_back(std::make_unique<GXFTrackState>(
4367 std::unique_ptr<const TrackParameters>(newparams)
4377 if (firsthit == firstmuonhit && cache.m_extmat && (firstcalopar !=
nullptr)) {
4378 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4380 if (cache.m_msEntrance ==
nullptr) {
4384 cache.m_msEntrance =
geometry->trackingVolume(
"MuonSpectrometerEntrance");
4390 if (cache.m_msEntrance ==
nullptr) {
4392 }
else if (cache.m_msEntrance->inside(firstcalopar->position())) {
4395 *cache.m_msEntrance,
4399 if (muonpar1 !=
nullptr) {
4407 rot.col(2) = trackdir;
4409 trans.linear().matrix() << rot;
4410 trans.translation() << muonpar1->
position() - .1 * trackdir;
4411 PlaneSurface curvlinsurf(trans);
4413 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4421 if (curvlinpar !=
nullptr) {
4422 muonpar1 = std::move(curvlinpar);
4426 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->clone());
4430 DistanceSolution distsol;
4432 if (muonpar1 !=
nullptr) {
4433 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4441 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4443 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4444 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4447 states[0]->associatedSurface(),
4451 matvec_used =
false;
4453 if (matvec && !matvec->empty()) {
4454 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4456 for (
int j = 0; j < (
int) matvec->size(); j++) {
4457 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4459 if (meb !=
nullptr) {
4463 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4464 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4467 !trajectory.m_straightline &&
4468 (meot->energyLoss() !=
nullptr) &&
4469 std::abs(meff->deltaE()) > 25 &&
4470 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4472 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4476 (*matvec)[j]->trackParameters() !=
nullptr
4477 ? (*matvec)[j]->trackParameters()->clone()
4480 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4482 std::unique_ptr<const TrackParameters>(tmpparams)
4495 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4496 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4499 newstates.reserve(oldstates.size() + matstates.size());
4502 int firstlayerno = -1;
4504 if (cache.m_acceleration) {
4505 trajectory.addBasicState(std::move(oldstates[0]));
4511 for (
int i = cache.m_acceleration ? 1 : 0;
i < (
int) oldstates.size();
i++) {
4512 bool addlayer =
true;
4514 while (addlayer && layerno < (
int) matstates.size()) {
4516 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4518 DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4519 layerpar->position(),
4520 layerpar->momentum().unit()
4525 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
4530 double cylinderradius = layerpar->associatedSurface().bounds().r();
4531 double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4533 if (trackimpact > cylinderradius - 5 *
mm) {
4539 if (
i == (
int) oldstates.size() - 1) {
4544 GXFMaterialEffects *meff = matstates[layerno]->materialEffects();
4546 if (meff->sigmaDeltaPhi() > .4 || (meff->sigmaDeltaPhi() == 0 && meff->sigmaDeltaE() <= 0)) {
4547 if (meff->sigmaDeltaPhi() > .4) {
4548 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4551 if (meff->sigmaDeltaPhi() == 0) {
4559 if (firstlayerno < 0) {
4560 firstlayerno = layerno;
4563 trajectory.addMaterialState(std::move(matstates[layerno]));
4565 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4566 const TrackingVolume *tvol =
m_navigator->volume(ctx,layerpar->position());
4567 const Layer *lay =
nullptr;
4569 if (tvol !=
nullptr) {
4570 lay = (tvol->closestMaterialLayer(layerpar->position(),layerpar->momentum().normalized())).
object;
4573 const MaterialProperties *matprop = lay !=
nullptr ? lay->fullUpdateMaterialProperties(*layerpar) :
nullptr;
4574 meff->setMaterialProperties(matprop);
4580 trajectory.addBasicState(std::move(oldstates[
i]));
4583 ATH_MSG_DEBUG(
"Total X0: " << trajectory.totalX0() <<
" total eloss: " << trajectory.totalEnergyLoss());
4585 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 3059 of file GlobalChi2Fitter.cxx.
3073 const double *
pos = parforextrap.position().data();
3075 cache.m_field_cache.getFieldZR(
pos,
field);
3081 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3082 double xc = parforextrap.position().x() -
r * sinphi;
3083 double yc = parforextrap.position().y() +
r * cosphi;
3084 double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
3085 double z0 = parforextrap.position().z();
3086 double d = xc * xc + yc * yc;
3087 double rcyl = surf.bounds().r();
3088 double mysqrt = ((
r + rcyl) * (
r + rcyl) -
d) * (
d - (
r - rcyl) * (
r - rcyl));
3094 mysqrt = std::sqrt(mysqrt);
3095 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3096 double secondterm = (mysqrt * yc) / (2 *
d);
3097 double x1 = firstterm + secondterm;
3098 double x2 = firstterm - secondterm;
3099 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3100 secondterm = (mysqrt * xc) / (2 *
d);
3101 double y1 = firstterm - secondterm;
3102 double y2 = firstterm + secondterm;
3103 double x = parforextrap.position().x();
3104 double y = parforextrap.position().y();
3105 double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3106 double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3108 if (dist1 < dist2) {
3116 double phi1 = std::atan2(
y - yc,
x - xc);
3117 double deltaphi = phi1 -
phi0;
3119 coercePhiCoordinateRange(deltaphi);
3121 double delta_z =
r * deltaphi / tantheta;
3122 double z =
z0 + delta_z;
3126 if (std::abs(
z - surf.center().z()) > surf.bounds().halflengthZ()) {
3131 double phidir = parforextrap.parameters()[
Trk::phi] + deltaphi;
3132 coercePhiCoordinateRange(phidir);
3136 double costracksurf = std::abs(normal.unit().dot(trackdir));
3138 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 3014 of file GlobalChi2Fitter.cxx.
3026 const double *
pos = parforextrap.position().data();
3028 cache.m_field_cache.getFieldZR(
pos,
field);
3036 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3037 double xc = parforextrap.position().x() -
r * sinphi;
3038 double yc = parforextrap.position().y() +
r * cosphi;
3039 double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
3040 double z0 = parforextrap.position().z();
3041 double delta_s = (surf.center().z() -
z0) / costheta;
3047 const DiscBounds *discbounds = (
const DiscBounds *) (&surf.bounds());
3049 if (
perp > discbounds->rMax() ||
perp < discbounds->rMin()) {
3053 double costracksurf = std::abs(costheta);
3055 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 3352 of file GlobalChi2Fitter.cxx.
3365 upstreamlayers.reserve(5);
3372 double firstz = firstsistate.trackParameters()->position().z();
3373 double firstr = firstsistate.trackParameters()->position().perp();
3374 double firstz2 = hasmat ? lastsistate.trackParameters()->position().z() : firstsistate.trackParameters()->position().z();
3375 double firstr2 = hasmat ? lastsistate.trackParameters()->position().perp() : firstsistate.trackParameters()->position().perp();
3377 GXFTrackState *firststate = oldstates.front().get();
3378 GXFTrackState *laststate = oldstates.back().get();
3384 double lastz = laststate->position().z();
3385 double lastr = laststate->position().perp();
3387 const Layer *startlayer = firststate->associatedSurface().associatedLayer();
3388 const Layer *startlayer2 = hasmat ? lastsistate.associatedSurface().associatedLayer() :
nullptr;
3389 const Layer *endlayer = laststate->associatedSurface().associatedLayer();
3392 double slope = (tantheta != 0) ? 1 / tantheta : 0;
3398 std::vector < const Layer *>::const_iterator
it;
3399 std::vector < const Layer *>::const_iterator itend;
3407 it = cache.m_posdiscs.begin();
3408 itend = cache.m_posdiscs.end();
3410 it = cache.m_negdiscs.begin();
3411 itend = cache.m_negdiscs.end();
3417 for (;
it != itend; ++
it) {
3422 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3430 const DiscBounds *discbounds = (
const DiscBounds *) (&(*it)->surfaceRepresentation().bounds());
3435 if (discbounds->rMax() < firstr || discbounds->rMin() > lastr) {
3439 double rintersect = firstr + ((*it)->surfaceRepresentation().center().z() - firstz) / slope;
3442 rintersect < discbounds->rMin() - 50 ||
3443 rintersect > discbounds->rMax() + 50
3453 if ((*
it) == endlayer) {
3465 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3468 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3475 (*
it) != startlayer &&
3476 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3477 (*it) == startlayer2)
3487 for (
const auto *barrelcylinder : cache.m_barrelcylinders) {
3491 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3498 double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().r() - firstr) * slope;
3500 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3501 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3505 if (barrelcylinder == endlayer) {
3512 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3513 barrelcylinder == startlayer) {
3514 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3517 if (barrelcylinder != startlayer &&
3518 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3519 barrelcylinder == startlayer2)) {
3520 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 3141 of file GlobalChi2Fitter.cxx.
3150 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3151 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(
states);
3162 for (
int i = 0;
i <= indexoffset;
i++) {
3163 trajectory.addBasicState(std::move(oldstates[
i]));
3172 for (
int i = indexoffset + 1;
i < (
int) oldstates.size();
i++) {
3173 double rmeas = oldstates[
i]->position().perp();
3174 double zmeas = oldstates[
i]->position().z();
3182 while (layerindex < (
int)
layers.size()) {
3184 double costracksurf = 0.0;
3199 const CylinderSurface *cylsurf = (
const CylinderSurface *) (&
layer->surfaceRepresentation());
3205 if (oldstates[
i]->trackParameters() !=
nullptr) {
3206 double rlayer = cylsurf->bounds().r();
3207 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->position().perp() - rlayer)) {
3208 parforextrap = oldstates[
i]->trackParameters();
3224 if (cylsurf->bounds().r() > rmeas)
break;
3232 const DiscSurface *discsurf = (
const DiscSurface *) (&
layer->surfaceRepresentation());
3234 if (oldstates[
i]->trackParameters() !=
nullptr) {
3235 double zlayer = discsurf->center().z();
3236 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->position().z() - zlayer)) {
3237 parforextrap = oldstates[
i]->trackParameters();
3248 if (std::abs(discsurf->center().z()) > std::abs(zmeas))
break;
3250 throw std::logic_error(
"Unhandled surface.");
3257 const MaterialProperties *matprop =
layer->layerMaterialProperties()->fullMaterial(
intersect);
3258 if (matprop ==
nullptr) {
3267 double X0 = matprop->thicknessInX0();
3271 double actualx0 =
X0 / costracksurf;
3272 double de = -std::abs(
3273 (matprop->thickness() / costracksurf) *
3276 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3279 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3281 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3285 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3286 meff->setDeltaE(de);
3287 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3288 meff->setX0(actualx0);
3289 meff->setSurface(&
layer->surfaceRepresentation());
3290 meff->setMaterialProperties(matprop);
3296 std::unique_ptr<EnergyLoss> eloss;
3298 if (cache.m_fiteloss || (matEffects ==
electron && cache.m_asymeloss)) {
3299 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3301 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3306 if (eloss !=
nullptr) {
3307 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3311 if (matEffects ==
electron && cache.m_asymeloss) {
3312 meff->setDeltaE(-5);
3313 if (trajectory.numberOfTRTHits() == 0) {
3314 meff->setScatteringSigmas(0, 0);
3317 meff->setSigmaDeltaE(50);
3318 if (eloss !=
nullptr) {
3319 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3320 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3325 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3326 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3327 " sigma eloss: " << meff->sigmaDeltaE()
3334 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3336 std::unique_ptr<const TrackParameters>()
3339 trajectory.addMaterialState(std::move(matstate));
3348 trajectory.addBasicState(std::move(oldstates[
i]));
◆ alignmentFit()
Definition at line 1862 of file GlobalChi2Fitter.cxx.
1868 const EventContext& ctx = Gaudi::Hive::currentContext();
1872 delete alignCache.m_derivMatrix;
1873 alignCache.m_derivMatrix =
nullptr;
1875 delete alignCache.m_fullCovarianceMatrix;
1876 alignCache.m_fullCovarianceMatrix =
nullptr;
1877 alignCache.m_iterationsOfLastFit = 0;
1880 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1881 if(newTrack !=
nullptr){
1882 if(cache.m_derivmat.size() != 0)
1883 alignCache.m_derivMatrix =
new Amg::MatrixX(cache.m_derivmat);
1884 if(cache.m_fullcovmat.size() != 0)
1885 alignCache.m_fullCovarianceMatrix =
new Amg::MatrixX(cache.m_fullcovmat);
1886 alignCache.m_iterationsOfLastFit = cache.m_lastiter;
◆ backupCombinationStrategy()
Definition at line 1279 of file GlobalChi2Fitter.cxx.
1287 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1289 bool firstismuon =
false;
1290 const Track *indettrack = &intrk1;
1293 indettrack = &intrk2;
1304 const auto *
const pParametersVector = indettrack->trackParameters();
1307 if (pParametersVector->size() > 1)
1308 firstidpar = (*pParametersVector)[1];
1310 firstidpar = pParametersVector->back();
1312 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1313 if ((firstidpar !=
nullptr) && (cache.m_caloEntrance !=
nullptr))
1317 if (lastidpar ==
nullptr) {
1318 lastidpar.reset(pParametersVector->back()->clone());
1321 std::unique_ptr < const TrackParameters > firstscatpar;
1322 std::unique_ptr < const TrackParameters > lastscatpar;
1323 std::unique_ptr < const TrackParameters > elosspar;
1328 lastidpar->position(),
1329 lastidpar->momentum(),
1331 lastidpar->position()
1338 calomeots[0].associatedSurface(),
1341 trajectory.m_fieldprop,
1344 if (!firstscatpar) {
1348 std::unique_ptr<const TrackParameters> tmppar(
1352 calomeots[1].associatedSurface(),
1355 trajectory.m_fieldprop,
1365 double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1366 double de = std::abs(calomeots[1].energyLoss()->deltaE());
1368 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1374 double newqoverp =
sign / std::sqrt(newp2);
1379 tmppar->associatedSurface().createUniqueTrackParameters(
1386 calomeots[2].associatedSurface(),
1389 trajectory.m_fieldprop,
1399 calomeots[2].associatedSurface(),
1402 trajectory.m_fieldprop,
1413 calomeots[1].associatedSurface(),
1416 trajectory.m_fieldprop,
1425 double newqoverp =
sign /
1426 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1427 std::abs(calomeots[1].energyLoss()->deltaE()));
1431 std::unique_ptr<const TrackParameters>tmppar(
1432 elosspar->associatedSurface().createUniqueTrackParameters(
1440 calomeots[0].associatedSurface(),
1443 trajectory.m_fieldprop,
1447 if (!firstscatpar) {
1452 for (; itStates != endState; ++itStates) {
1460 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1462 cache.m_idmat =
false;
1473 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1474 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1475 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1479 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1480 elossmeff->setSigmaDeltaE(sigmadp);
1483 elossmeff->setdelta_p(
dp);
1485 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1486 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1487 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1489 GXFTrackState *secondscatstate = trajectory.trackStates().back().get();
1490 const Surface *triggersurf1 =
nullptr;
1491 const Surface *triggersurf2 =
nullptr;
1495 bool seenmdt =
false;
1496 bool mdtbetweenphihits =
false;
1500 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1501 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1502 (!firstismuon ? ++itStates2 : --itStates2)
1505 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1510 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1511 const Surface *surf = &pMeasurement->associatedSurface();
1515 if (isCompetingRIOsOnTrack) {
1517 rot = &crot->rioOnTrack(0);
1520 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1523 if ((rot !=
nullptr) &&
m_DetID->
is_mdt(rot->identify()) && (triggersurf1 !=
nullptr)) {
1527 (rot !=
nullptr) && (
1533 bool measphi =
true;
1534 Amg::Vector3D measdir = surf->transform().rotation().col(0);
1536 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1537 if (std::abs(dotprod1) > .5 || std::abs(dotprod2) > .5) {
1543 (*itStates2)->trackParameters() !=
nullptr ?
1544 (*itStates2)->trackParameters()->position() :
1545 rot->globalPosition();
1546 if (triggersurf1 !=
nullptr) {
1547 triggerpos2 = thispos;
1548 triggersurf2 = surf;
1550 mdtbetweenphihits =
true;
1553 triggerpos1 = thispos;
1554 triggersurf1 = surf;
1560 double mdttrig1 = 999999;
1561 double mdttrig2 = 999999;
1562 const Surface *mdtsurf1 =
nullptr;
1563 const Surface *mdtsurf2 =
nullptr;
1566 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1567 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1568 (!firstismuon ? ++itStates2 : --itStates2)
1570 const Surface *surf =
nullptr;
1572 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1575 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1578 if (surf ==
nullptr) {
1581 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1585 if (isCompetingRioOnTrack) {
1587 rot = &crot->rioOnTrack(0);
1590 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1593 const bool thisismdt = rot and
m_DetID->
is_mdt(rot->identify());
1596 (*itStates2)->trackParameters() !=
nullptr ?
1597 (*itStates2)->trackParameters()->position() :
1598 pThisMeasurement->globalPosition();
1599 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).
mag() < mdttrig1) {
1600 mdttrig1 = (globpos - triggerpos1).
mag();
1603 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).
mag() < mdttrig2) {
1604 mdttrig2 = (globpos - triggerpos2).
mag();
1610 GXFTrackState * firstpseudostate =
nullptr;
1611 std::vector<GXFTrackState *> outlierstates;
1612 std::vector<GXFTrackState *> outlierstates2;
1614 outlierstates.reserve(10);
1615 outlierstates2.reserve(10);
1617 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1619 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1620 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1622 bool isStraightLine =
1623 pMeasurement !=
nullptr ?
1630 (newpseudo ==
nullptr) && (
1631 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1632 std::abs(pMeasurement->globalPosition().z()) < 10000
1635 std::unique_ptr<const TrackParameters> par2;
1636 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1637 par2.reset((*itStates2)->trackParameters()->clone());
1641 *secondscatstate->trackParameters(),
1642 pMeasurement->associatedSurface(),
1645 trajectory.m_fieldprop,
1649 if (par2 ==
nullptr) {
1655 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1658 par2->associatedSurface()
1661 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1668 firstpseudo->setMeasurementErrors(
errors);
1669 firstpseudostate = firstpseudo.get();
1670 trajectory.addMeasurementState(std::move(firstpseudo));
1675 if (isPseudo && !firstismuon) {
1679 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1681 cache.m_idmat =
false;
1689 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1690 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1691 (mdtsurf1 !=
nullptr)
1693 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->transform());
1695 transf->translation() << triggerpos1;
1696 StraightLineSurface slsurf(*transf);
1700 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1704 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1711 pseudostate1->setMeasurementErrors(
errors);
1712 outlierstates2.push_back(pseudostate1.get());
1713 trajectory.addMeasurementState(std::move(pseudostate1));
1717 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1718 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1719 mdtbetweenphihits &&
1720 (mdtsurf2 !=
nullptr)
1722 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->transform());
1723 transf->translation() << triggerpos2;
1724 StraightLineSurface slsurf(*transf);
1728 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1732 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1739 pseudostate2->setMeasurementErrors(
errors);
1741 outlierstates2.push_back(pseudostate2.get());
1742 trajectory.addMeasurementState(std::move(pseudostate2));
1749 trajectory.trackStates().back()->measurementType() ==
TrackState::TGC ||
1751 trajectory.trackStates().back()->measurementType() ==
TrackState::RPC &&
1752 trajectory.trackStates().back()->measuresPhi()
1757 outlierstates.push_back(trajectory.trackStates().back().get());
1758 trajectory.setOutlier((
int) trajectory.trackStates().size() - 1,
true);
1763 trajectory.setNumberOfPerigeeParameters(0);
1767 trajectory.setPrefit(2);
1769 cache.m_matfilled =
true;
1770 bool tmpacc = cache.m_acceleration;
1771 cache.m_acceleration =
false;
1773 std::unique_ptr<Trk::Track> tmp_track(
1774 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1775 cache.m_acceleration = tmpacc;
1777 cache.m_matfilled =
false;
1780 trajectory.converged() &&
1781 std::abs(trajectory.residuals().tail<1>()(0) / trajectory.errors().tail<1>()(0)) > 10
1786 if (trajectory.converged()) {
1787 if (firstpseudostate !=
nullptr) {
1792 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1795 par2->associatedSurface()
1797 firstpseudostate->setMeasurement(std::move(newpseudo));
1798 firstpseudostate->setRecalibrated(
false);
1801 for (
int j = 0; j < (
int) trajectory.trackStates().size(); j++) {
1802 for (
const auto &
i : outlierstates2) {
1803 if (trajectory.trackStates()[j].get() ==
i) {
1804 trajectory.setOutlier(j,
true);
1808 for (
const auto &
i : outlierstates) {
1809 if (trajectory.trackStates()[j].get() ==
i) {
1810 trajectory.setOutlier(j,
false);
1816 itStates = (firstismuon ? beginStates2 : endState - 1);
1817 itStates != (firstismuon ? endState2 : beginStates - 1);
1818 (firstismuon ? ++itStates : --itStates)
1824 makeProtoState(cache, trajectory, *itStates, (firstismuon ? -1 : 0));
1828 trajectory.setPrefit(0);
1829 trajectory.setNumberOfPerigeeParameters(5);
1831 cache.m_matfilled =
false;
◆ calculateDerivatives()
void Trk::GlobalChi2Fitter::calculateDerivatives |
( |
GXFTrajectory & |
trajectory | ) |
|
|
staticprivate |
Definition at line 7817 of file GlobalChi2Fitter.cxx.
7818 int nstatesupstream = trajectory.numberOfUpstreamStates();
7819 int nscatupstream = trajectory.numberOfUpstreamScatterers();
7820 int nbremupstream = trajectory.numberOfUpstreamBrems();
7821 int nscats = trajectory.numberOfScatterers();
7822 int nperpars = trajectory.numberOfPerigeeParameters();
7823 int nfitpars = trajectory.numberOfFitParameters();
7825 using Matrix55 = Eigen::Matrix<double, 5, 5>;
7827 Matrix55 initialjac;
7828 initialjac.setZero();
7829 initialjac(4, 4) = 1;
7831 Matrix55 jacvertex(initialjac);
7833 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.numberOfScatterers(), initialjac);
7834 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.numberOfBrems(), initialjac);
7836 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7837 GXFTrackState *prevstate =
nullptr, *state =
nullptr;
7839 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
7841 for (
bool forward : {
false,
true}) {
7843 hit_begin = nstatesupstream;
7845 scatno = nscatupstream;
7846 bremno = nbremupstream;
7848 hit_begin = nstatesupstream - 1;
7850 scatno = trajectory.numberOfUpstreamScatterers() - 1;
7851 bremno = trajectory.numberOfUpstreamBrems() - 1;
7855 int hitno = hit_begin;
7856 forward ? (hitno < hit_end) : (hitno >= hit_end);
7857 hitno += (forward ? 1 : -1)
7860 state =
states[hitno].get();
7864 if (fillderivmat && state->derivatives().cols() != nfitpars) {
7865 state->derivatives().resize(5, nfitpars);
7866 state->derivatives().setZero();
7869 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
7871 if (hitno == (forward ? hit_end - 1 : 0)) {
7872 if (!fillderivmat) {
7880 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
7882 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
7883 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
7884 jacvertex(4, 4) = jac(4, 4);
7886 int jmin = 0, jmax = 0, jcnt = 0;
7887 int lp_bgn = 0, lp_end = 0;
7891 jcnt = jmax - jmin + 1;
7893 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
7896 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7898 i == scatno + (forward ? -1 : 1) &&
7899 prevstate !=
nullptr &&
7901 (!trajectory.prefit() || prevstate->materialEffects()->deltaE() == 0)
7903 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7904 jacscat[
i](4, 4) = jac(4, 4);
7906 calculateJac(jac, jacscat[
i], jmin, jmax);
7910 Eigen::MatrixXd & derivmat = state->derivatives();
7911 int scatterPos = nperpars + 2 *
i;
7913 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
7919 jcnt = jmax - jmin + 1;
7921 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
7924 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7926 i == bremno + (forward ? -1 : 1) &&
7928 prevstate->materialEffects() &&
7929 prevstate->materialEffects()->sigmaDeltaE() > 0
7931 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7932 jacbrem[
i](4, 4) = jac(4, 4);
7934 calculateJac(jac, jacbrem[
i], jmin, jmax);
7938 Eigen::MatrixXd & derivmat = state->derivatives();
7939 int scatterPos = nperpars + 2 * nscats +
i;
7941 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
7945 calculateJac(jac, jacvertex, 0, 4);
7949 Eigen::MatrixXd & derivmat = state->derivatives();
7950 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
7952 if (nperpars == 5) {
7953 derivmat.col(4).segment(0, 4) *= .001;
7954 derivmat(4, 4) = .001 * jacvertex(4, 4);
7960 (!trajectory.prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
7962 scatno += (forward ? 1 : -1);
7966 states[hitno]->materialEffects() &&
7967 states[hitno]->materialEffects()->sigmaDeltaE() > 0
7969 bremno += (forward ? 1 : -1);
7972 prevstate =
states[hitno].get();
◆ calculateTrackErrors()
Definition at line 7979 of file GlobalChi2Fitter.cxx.
7987 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7988 int nstatesupstream = trajectory.numberOfUpstreamStates();
7990 GXFTrackState *prevstate =
nullptr;
7991 int i = nstatesupstream;
7992 for (
int j = 0; j < (
int)
states.size(); j++) {
7993 if (j < nstatesupstream) {
8000 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8001 if (stateno == 0 || stateno == nstatesupstream) {
8002 prevstate =
nullptr;
8005 std::unique_ptr<GXFTrackState> & state =
states[
index];
8006 if (state->materialEffects() !=
nullptr) {
8007 prevstate = state.get();
8011 if (!state->hasTrackCovariance()) {
8012 state->zeroTrackCovariance();
8014 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8016 if ((prevstate !=
nullptr) &&
8020 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8023 trackerrmat = jac * prevcov * jac.transpose();
8027 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8031 const MeasurementBase *measurement = state->measurement();
8032 const Amg::MatrixX & meascov = measurement->localCovariance();
8034 ParamDefsAccessor paraccessor;
8038 bool errorok =
true;
8039 for (
int i = 0;
i < 5;
i++) {
8040 if (measurement->localParameters().contains(paraccessor.pardef[
i])) {
8042 && trackerrmat(
i,
i) > meascov(j, j)) {
8044 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8045 trackerrmat(
i,
i) = meascov(j, j);
8046 for (
int k = 0;
k < 5;
k++) {
8056 for (
int i = 0;
i < 5;
i++) {
8060 for (
int j = 0; j < 5; j++) {
8067 if (trajectory.m_straightline) {
8068 trackerrmat(4, 4) = 1
e-20;
8072 state->trackParameters();
8074 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8076 if (state->hasTrackCovariance()) {
8077 trkerrmat = (state->trackCovariance());
8079 trkerrmat = std::nullopt;
8082 const AmgVector(5) & tpars = tmptrackpar->parameters();
8083 std::unique_ptr<const TrackParameters> trackpar(
8084 tmptrackpar->associatedSurface().createUniqueTrackParameters(tpars[0],
8089 std::move(trkerrmat))
8091 state->setTrackParameters(std::move(trackpar));
8092 FitQualityOnSurface fitQual{};
8094 if (errorok && trajectory.nDOF() > 0) {
8095 fitQual =
m_updator->fullStateFitQuality(
8096 *state->trackParameters(),
8097 measurement->localParameters(),
8098 measurement->localCovariance()
8101 fitQual = FitQualityOnSurface(0, state->numberOfMeasuredParameters());
8104 state->setFitQuality(fitQual);
8106 prevstate = state.get();
◆ calculateTrackParameters()
Definition at line 7593 of file GlobalChi2Fitter.cxx.
7601 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7602 int nstatesupstream = trajectory.numberOfUpstreamStates();
7603 const TrackParameters *prevtrackpar = trajectory.referenceParameters();
7604 std::unique_ptr<const TrackParameters> tmptrackpar;
7606 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7610 DistanceSolution distsol = surf1.straightLineDistanceEstimate(
7611 prevtrackpar->position(), prevtrackpar->momentum().unit()
7618 distsol.numberOfSolutions() > 0 &&
7619 prevtrackpar != trajectory.referenceParameters()
7629 trajectory.m_fieldprop,
7636 (
rv.m_parameters !=
nullptr) &&
7637 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7643 if (
rv.m_parameters ==
nullptr) {
7644 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7645 " pos: " << prevtrackpar->position() <<
" destination surface: " << surf1);
7649 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7653 if (
rv.m_jacobian != std::nullopt) {
7655 states[hitno]->materialEffects() !=
nullptr &&
7656 states[hitno]->materialEffects()->deltaE() != 0 &&
7657 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7658 !trajectory.m_straightline
7660 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7661 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7662 double mass = trajectory.mass();
7663 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7664 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7667 states[hitno]->setJacobian(*
rv.m_jacobian);
7668 }
else if (calcderiv) {
7673 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7675 if (meff !=
nullptr && hitno != 0) {
7676 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7677 surf, *meff, *
states[hitno]->trackParameters(), trajectory.mass(), -1
7680 if (std::holds_alternative<FitterStatusCode>(
r)) {
7681 return std::get<FitterStatusCode>(
r);
7684 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7685 prevtrackpar = tmptrackpar.get();
7687 prevtrackpar = currenttrackpar;
7691 prevtrackpar = trajectory.referenceParameters();
7693 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7696 DistanceSolution distsol = surf.straightLineDistanceEstimate(prevtrackpar->position(), prevtrackpar->momentum().unit());
7700 if (
distance < 0 && distsol.numberOfSolutions() > 0 && prevtrackpar != trajectory.referenceParameters()) {
7709 trajectory.m_fieldprop,
7715 (
rv.m_parameters !=
nullptr) &&
7717 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7722 if (
rv.m_parameters ==
nullptr) {
7723 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7724 " pos: " << prevtrackpar->
7725 position() <<
" destination surface: " << surf);
7729 if (
rv.m_jacobian != std::nullopt) {
7731 states[hitno]->materialEffects() !=
nullptr &&
7732 states[hitno]->materialEffects()->deltaE() != 0 &&
7733 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7734 !trajectory.m_straightline
7736 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7737 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7738 double mass = trajectory.mass();
7739 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7742 newp = std::sqrt(newp);
7745 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7748 states[hitno]->setJacobian(*
rv.m_jacobian);
7749 }
else if (calcderiv) {
7754 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7756 if (meff !=
nullptr) {
7757 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7758 surf, *meff, *
rv.m_parameters, trajectory.mass(), +1
7761 if (std::holds_alternative<FitterStatusCode>(
r)) {
7762 return std::get<FitterStatusCode>(
r);
7765 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7768 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7769 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 7567 of file GlobalChi2Fitter.cxx.
7576 PropagationResult
rv;
7579 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7582 if (
rv.m_parameters ==
nullptr) {
7583 propdir = invertPropdir(propdir);
7586 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 7528 of file GlobalChi2Fitter.cxx.
7537 std::unique_ptr<const TrackParameters>
rv;
7538 std::optional<TransportJacobian> jac{};
7549 if (
rv !=
nullptr && calcderiv) {
7554 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7560 return PropagationResult {
7563 std::move(extrapolation)
◆ correctAngles()
bool Trk::GlobalChi2Fitter::correctAngles |
( |
double & |
phi, |
|
|
double & |
theta |
|
) |
| |
|
staticprivate |
◆ fillDerivatives()
void Trk::GlobalChi2Fitter::fillDerivatives |
( |
GXFTrajectory & |
traj, |
|
|
bool |
onlybrem = false |
|
) |
| const |
|
private |
Definition at line 5498 of file GlobalChi2Fitter.cxx.
5504 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5508 int nscatupstream = trajectory.numberOfUpstreamScatterers();
5509 int nbremupstream = trajectory.numberOfUpstreamBrems();
5510 int nscat = trajectory.numberOfScatterers();
5511 int nbrem = trajectory.numberOfBrems();
5512 int nperparams = trajectory.numberOfPerigeeParameters();
5514 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5517 int nmeas = (
int) weightderiv.rows();
5519 ParamDefsAccessor paraccessor;
5521 for (std::unique_ptr<GXFTrackState> & state :
states) {
5524 ((state->materialEffects() ==
nullptr) || state->materialEffects()->sigmaDeltaE() <= 0)
5530 const MeasurementBase *measbase = state->measurement();
5531 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5532 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5535 double sinstereo = 0;
5538 sinstereo = state->sinStereo();
5541 double cosstereo = (sinstereo == 0) ? 1. : std::sqrt(1 - sinstereo * sinstereo);
5543 for (
int i = 0;
i < 5;
i++) {
5545 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5551 if (trajectory.numberOfPerigeeParameters() > 0) {
5552 int cols = trajectory.m_straightline ? 4 : 5;
5555 weightderiv.row(measno).head(
cols) =
5556 (derivatives.row(0).head(
cols) * cosstereo +
5557 sinstereo * derivatives.row(1).head(
cols)) /
5560 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5564 for (
int j = scatmin; j < scatmax; j++) {
5565 int index = nperparams + ((trajectory.prefit() != 1) ? 2 * j : j);
5566 double thisderiv = 0;
5570 if (
i == 0 && sinstereo != 0) {
5571 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5576 weightderiv(measno,
index) = thisderiv /
error[measno];
5578 if (trajectory.prefit() != 1) {
5581 if (
i == 0 && sinstereo != 0) {
5582 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5587 weightderiv(measno,
index) = thisderiv /
error[measno];
5591 for (
int j = bremmin; j < bremmax; j++) {
5592 double thisderiv = 0;
5593 int index = j + nperparams + 2 * nscat;
5594 if (
i == 0 && sinstereo != 0) {
5595 thisderiv = derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index);
5597 thisderiv = derivatives(
i,
index);
5599 weightderiv(measno,
index) = thisderiv /
error[measno];
5605 double *
errors = state->measurementErrors();
5606 for (
int i = 0;
i < 5;
i++) {
5613 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5618 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5620 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5621 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5623 double mass = .001 * trajectory.mass();
5625 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5629 auto multiplier = [] (
double mass,
double qOverP){
5632 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5633 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5635 if (trajectory.numberOfPerigeeParameters() > 0) {
5636 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5639 const auto bremNoBase = nperparams + 2 * nscat;
5640 if (bremno < nbremupstream) {
5641 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5642 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5643 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5646 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5647 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5648 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
◆ fillResiduals()
Definition at line 5188 of file GlobalChi2Fitter.cxx.
5200 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5205 int nbrem = trajectory.numberOfBrems();
5206 int nperpars = trajectory.numberOfPerigeeParameters();
5207 int nfitpars = trajectory.numberOfFitParameters();
5210 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5211 int nidhits = trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits();
5212 int nsihits = trajectory.numberOfSiliconHits();
5213 int ntrthits = trajectory.numberOfTRTHits();
5214 int nhits = trajectory.numberOfHits();
5215 int nmeas = (
int)
res.size();
5217 ParamDefsAccessor paraccessor;
5218 bool scatwasupdated =
false;
5220 GXFTrackState *state_maxbrempull =
nullptr;
5221 int bremno_maxbrempull = 0;
5222 double maxbrempull = 0;
5224 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5225 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5228 const MeasurementBase *measbase = state->measurement();
5235 trajectory.nDOF() != 0 &&
5236 std::abs((trajectory.prevchi2() - trajectory.chi2()) / trajectory.nDOF()) < 15 &&
5237 !state->associatedSurface().isFree() &&
5238 nidhits < trajectory.numberOfHits() &&
5239 (nperpars == 0 || nidhits > 0) &&
5240 (!state->isRecalibrated())
5245 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5248 currenttrackpar->associatedSurface()
5251 state->setMeasurement(std::move(newpseudo));
5252 measbase = state->measurement();
5255 double *
errors = state->measurementErrors();
5259 for (
int i = 0;
i < 5;
i++) {
5261 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5272 res[measno] = residuals[
i];
5274 if (
i == 2 && std::abs(std::abs(
res[measno]) - 2 *
M_PI) < std::abs(
res[measno])) {
5275 if (
res[measno] < 0) {
5284 double *
errors = state->measurementErrors();
5285 for (
int i = 0;
i < 5;
i++) {
5295 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5297 double deltaphi = state->materialEffects()->deltaPhi();
5298 double measdeltaphi = state->materialEffects()->measuredDeltaPhi();
5299 double sigmadeltaphi = state->materialEffects()->sigmaDeltaPhi();
5300 double deltatheta = state->materialEffects()->deltaTheta();
5301 double sigmadeltatheta = state->materialEffects()->sigmaDeltaTheta();
5303 if (trajectory.prefit() != 1) {
5304 b[nperpars + 2 * scatno] -= (deltaphi - measdeltaphi) / (sigmadeltaphi * sigmadeltaphi);
5305 b[nperpars + 2 * scatno + 1] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5307 b[nperpars + scatno] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5311 deltaphi * deltaphi / (sigmadeltaphi * sigmadeltaphi) +
5312 deltatheta * deltatheta / (sigmadeltatheta * sigmadeltatheta)
5318 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5319 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5321 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5322 const double pbrem = 1. / std::abs(qoverpbrem);
5323 const double p = 1. / std::abs(qoverp);
5324 const double mass = .001 * trajectory.mass();
5326 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5328 res[nmeas - nbrem + bremno] = .001 * averagenergyloss -
energy + bremEnergy;
5330 double sigde = state->materialEffects()->sigmaDeltaE();
5331 double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5332 double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5334 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5336 if (state->materialEffects()->isKink()) {
5337 maxbrempull = -999999999;
5338 state_maxbrempull =
nullptr;
5342 cache.m_asymeloss &&
5344 (trajectory.prefit() == 0) &&
5346 sigde != sigdepos &&
5349 double elosspull =
res[nmeas - nbrem + bremno] / (.001 * sigde);
5351 if (trajectory.mass() > 100) {
5352 if (elosspull < -1) {
5353 state->materialEffects()->setSigmaDeltaE(sigdepos);
5354 }
else if (elosspull > 1) {
5355 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5358 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5359 }
else if ((trajectory.numberOfTRTHits() == 0) ||
it >= 3) {
5361 !state->materialEffects()->isKink() && (
5362 (elosspull < -.2 &&
m_fixbrem == -1 && elosspull < maxbrempull) ||
5366 bremno_maxbrempull = bremno;
5367 state_maxbrempull = state.get();
5368 maxbrempull = elosspull;
5377 (trajectory.prefit() == 0) &&
5378 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5379 state->materialEffects()->isMeasuredEloss() &&
5380 res[nmeas - nbrem + bremno] / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5383 trajectory.prefit() != 0 ? trajectory.referenceParameters()
5384 :
states[hitno - 2]->trackParameters();
5387 std::vector<MaterialEffectsOnTrack> calomeots =
5392 parforcalo->associatedSurface(),
5396 if (calomeots.size() == 3) {
5397 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5398 double newres = .001 * averagenergyloss -
energy + bremEnergy;
5399 double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5401 if (std::abs(newres / newerr) < std::abs(
res[nmeas - nbrem + bremno] /
error[nmeas - nbrem + bremno])) {
5402 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5404 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5405 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5406 res[nmeas - nbrem + bremno] = newres;
5407 error[nmeas - nbrem + bremno] = newerr;
5411 state->materialEffects()->setMeasuredEloss(
false);
5419 for (; measno < nmeas; measno++) {
5420 if (
error[measno] == 0) {
5427 if (!doderiv && (scatwasupdated)) {
5431 double oldchi2 = trajectory.chi2();
5432 trajectory.setPrevChi2(oldchi2);
5433 trajectory.setChi2(
chi2);
5435 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5436 double newredchi2 = (trajectory.nDOF() > 0) ?
chi2 / trajectory.nDOF() : 0;
5439 trajectory.prefit() > 0 && (
5440 (newredchi2 < 2 &&
it != 0) ||
5441 (newredchi2 < oldredchi2 + .1 && std::abs(newredchi2 - oldredchi2) < 1 &&
it != 1)
5444 trajectory.setConverged(
true);
5447 double maxdiff = (nsihits != 0 && nsihits + ntrthits == nhits &&
chi2 < oldchi2) ? 200 : 1.;
5449 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5451 if (miniter < cache.m_miniter) {
5452 miniter = cache.m_miniter;
5456 trajectory.setConverged(
true);
5459 if ((state_maxbrempull !=
nullptr) && trajectory.converged()) {
5460 state_maxbrempull->materialEffects()->setSigmaDeltaE(
5461 10 * state_maxbrempull->materialEffects()->sigmaDeltaEPos()
5464 state_maxbrempull->materialEffects()->setKink(
true);
5465 trajectory.setConverged(
false);
5467 double olderror =
error[nmeas - nbrem + bremno_maxbrempull];
5468 double newerror = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5469 error[nmeas - nbrem + bremno_maxbrempull] = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5471 if (
a.cols() != nfitpars) {
5475 for (
int i = 0;
i < nfitpars;
i++) {
5476 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5480 for (
int j =
i; j < nfitpars; j++) {
5484 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5485 weightderiv(nmeas - nbrem + bremno_maxbrempull, j) *
5486 (1 - olderror * olderror / (newerror * newerror))
5490 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= olderror / newerror;
5493 trajectory.setChi2(1e15);
◆ finalize()
StatusCode Trk::GlobalChi2Fitter::finalize |
( |
| ) |
|
|
overridevirtual |
Definition at line 286 of file GlobalChi2Fitter.cxx.
289 if (m_fit_status[
S_FITS] > 0) {
292 <<
" track fits failed because of a matrix inversion failure");
294 <<
" tracks were rejected by the outlier logic");
296 <<
" track fits failed because of a propagation failure");
298 <<
" track fits failed because of an invalid angle (theta/phi)");
300 <<
" track fits failed because the fit did not converge");
302 <<
" tracks did not pass the chi^2 cut");
304 <<
" tracks were killed by the energy loss update");
307 return StatusCode::SUCCESS;
◆ fit() [1/6]
Definition at line 2445 of file GlobalChi2Fitter.cxx.
2452 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2457 GXFTrajectory trajectory;
2460 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2465 for (
const auto *itSet : rots) {
2466 if (itSet ==
nullptr) {
2467 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2473 std::unique_ptr<const TrackParameters> startpar(param.clone());
2476 matEffects ==
muon &&
2477 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0
2479 cache.m_matfilled =
true;
2480 trajectory.setPrefit(2);
2482 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2484 cache.m_matfilled =
false;
2486 if (!trajectory.converged()) {
2490 trajectory.setConverged(
false);
2491 const TrackParameters *firstpar = trajectory.trackStates()[0]->trackParameters();
2492 const TrackParameters *lastpar = trajectory.trackStates().back()->trackParameters();
2494 PerigeeSurface persurf(firstpar->position() - 10 * firstpar->momentum().unit());
2500 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2503 firstpar->associatedSurface()
2506 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2513 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2516 lastpar->associatedSurface()
2519 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2522 if (!trajectory.m_straightline) {
2523 trajectory.setPrefit(3);
2524 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2525 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2526 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2531 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2533 cache.m_matfilled =
true;
2535 if (!trajectory.converged()) {
2540 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2541 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2542 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2546 trajectory.setPrefit(0);
2549 firstpar = trajectory.trackStates().front()->trackParameters();
2554 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2557 firstpar->associatedSurface()
2560 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2564 trajectory.trackStates().front()->setMeasurementErrors(
errors);
2568 lastpar = trajectory.trackStates().back()->trackParameters();
2573 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2576 lastpar->associatedSurface()
2579 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2583 trajectory.trackStates().back()->setMeasurementErrors(
errors);
2587 std::unique_ptr<Track>
track;
2589 if (startpar !=
nullptr) {
2590 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2593 if (
track !=
nullptr) {
2596 cache.m_matfilled =
false;
◆ fit() [2/6]
Definition at line 2223 of file GlobalChi2Fitter.cxx.
2230 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2233 for (
const auto *prd : prds) {
2234 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2236 const PlaneSurface *plsurf =
nullptr;
2239 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2241 const StraightLineSurface *slsurf =
nullptr;
2244 slsurf =
static_cast < const StraightLineSurface *
>(&prdsurf);
2246 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2247 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2252 }
else if (slsurf !=
nullptr) {
2261 }
else if (plsurf !=
nullptr) {
2262 if (param.covariance() !=
nullptr) {
2284 if (rot !=
nullptr) {
2285 rots.push_back(rot);
2289 std::unique_ptr<Track>
track =
2290 fit(ctx, rots, param, runOutlier, matEffects);
2292 for (
const auto *rot : rots) {
◆ fit() [3/6]
Definition at line 2300 of file GlobalChi2Fitter.cxx.
2307 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2312 GXFTrajectory trajectory;
2315 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2322 if (minpar ==
nullptr) {
2323 minpar = *(inputTrack.trackParameters()->begin());
2332 bool old_reintoutl = cache.m_reintoutl;
2333 cache.m_reintoutl =
false;
2334 bool tmpasymeloss = cache.m_asymeloss;
2337 cache.m_asymeloss =
true;
2340 for (; itStates != endState; ++itStates) {
2344 !trajectory.trackStates().empty() &&
2345 (trajectory.trackStates().back()->materialEffects() !=
nullptr) &&
2346 trajectory.trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2348 trajectory.trackStates().back()->materialEffects()->setKink(
true);
2352 cache.m_reintoutl = old_reintoutl;
2353 MeasurementSet::const_iterator itSet = addMeasColl.begin();
2354 MeasurementSet::const_iterator itSetEnd = addMeasColl.end();
2356 for (; itSet != itSetEnd; ++itSet) {
2357 if ((*itSet) ==
nullptr) {
2358 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2365 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2366 cache.m_asymeloss = tmpasymeloss;
2368 if (
track !=
nullptr) {
2370 inputTrack.fitQuality()->numberDoF() != 0 ?
2371 inputTrack.fitQuality()->chiSquared() / inputTrack.fitQuality()->numberDoF() :
2375 track->fitQuality()->numberDoF() != 0 ?
2376 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() :
2379 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2383 track.reset(
nullptr);
2387 if (
track !=
nullptr) {
◆ fit() [4/6]
Definition at line 2397 of file GlobalChi2Fitter.cxx.
2404 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2408 for (
const auto *prd : prds) {
2409 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2412 std::unique_ptr<const TrackParameters>trackparForCorrect(
2413 hitparam->associatedSurface().createUniqueTrackParameters(
2428 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect);
2431 if (rot !=
nullptr) {
2432 rots.push_back(rot);
2436 std::unique_ptr<Track>
track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2438 for (
const auto *rot : rots) {
◆ fit() [5/6]
Definition at line 1838 of file GlobalChi2Fitter.cxx.
1844 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1849 GXFTrajectory trajectory;
1852 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1857 return std::unique_ptr<Track>(
1858 fitIm(ctx, cache, inputTrack, runOutlier, matEffects));
◆ fit() [6/6]
Definition at line 313 of file GlobalChi2Fitter.cxx.
320 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
325 GXFTrajectory trajectory;
327 trajectory.m_straightline = (
328 !cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn()
335 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
336 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
338 bool measphi =
false;
340 for (
const auto *
i : *(muontrack->measurementsOnTrack())) {
345 rot = &crot->rioOnTrack(0);
352 const Surface *surf = &rot->associatedSurface();
355 double dotprod2 = measdir.dot(
357 surf->center().perp());
358 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
366 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
368 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
372 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
374 if (!cache.m_field_cache.solenoidOn()) {
375 const AmgVector(5) & newpars = parforcalo->parameters();
377 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
378 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
382 std::vector < MaterialEffectsOnTrack > calomeots;
385 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
389 parforcalo->associatedSurface(),
402 if (calomeots.empty()) {
407 std::unique_ptr<Track>
track;
410 cache.m_calomat =
false;
411 bool tmp2 = cache.m_extmat;
412 bool tmp4 = cache.m_idmat;
417 double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
418 double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
420 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
421 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
425 (errmatid !=
nullptr) &&
426 (errmatmuon !=
nullptr) &&
431 double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
432 double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
433 double energyerror = std::sqrt(
434 calomeots[1].energyLoss()->sigmaDeltaE() *
435 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
436 pmuonerror * pmuonerror
440 (std::abs(calomeots[1].energyLoss()->deltaE()) -
441 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
444 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
449 parforcalo->associatedSurface(),
454 if (calomeots.empty()) {
461 int nfits = cache.m_fit_status[
S_FITS];
462 bool firstfitwasattempted =
false;
464 if (cache.m_caloEntrance ==
nullptr) {
468 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
471 if (cache.m_caloEntrance ==
nullptr) {
478 (!cache.m_field_cache.toroidOn() && !cache.m_field_cache.solenoidOn()) ||
480 cache.m_getmaterialfromtrack &&
484 qoverpid * qoverpmuon > 0
489 if (cache.m_fit_status[
S_FITS] == (
unsigned int) (nfits + 1)) {
490 firstfitwasattempted =
true;
495 (
track ==
nullptr) &&
496 !firstfitwasattempted &&
497 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn())
500 GXFTrajectory trajectory2;
501 trajectory2.m_straightline = trajectory.m_straightline;
502 trajectory2.m_fieldprop = trajectory.m_fieldprop;
503 trajectory = trajectory2;
507 bool pseudoupdated =
false;
509 if (
track !=
nullptr) {
510 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
511 if (pseudostate ==
nullptr) {
522 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
527 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
529 pseudostate->measurement()->localParameters(),
530 pseudostate->measurement()->localCovariance()
533 if (updpar ==
nullptr) {
540 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
545 pseudopar->associatedSurface()
548 pseudostate->setMeasurement(std::move(newpseudo));
552 pseudostate->setMeasurementErrors(
errors);
553 pseudoupdated =
true;
557 trajectory.setConverged(
false);
558 cache.m_matfilled =
true;
564 *
track->perigeeParameters(),
566 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
569 cache.m_matfilled =
false;
573 cache.m_fit_status[
S_FITS] = nfits + 1;
575 if (
track !=
nullptr) {
576 track->info().addPatternReco(intrk1.info());
577 track->info().addPatternReco(intrk2.info());
581 cache.m_calomat =
tmp;
582 cache.m_extmat =
tmp2;
583 cache.m_idmat = tmp4;
◆ fitIm()
Definition at line 1892 of file GlobalChi2Fitter.cxx.
1900 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1902 GXFTrajectory trajectory;
1905 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1910 if (inputTrack.trackStateOnSurfaces()->empty()) {
1915 if (inputTrack.trackParameters()->empty()) {
1916 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1920 std::unique_ptr<const TrackParameters> minpar =
unique_clone(inputTrack.perigeeParameters());
1924 if (minpar ==
nullptr) {
1925 minpar =
unique_clone(*(inputTrack.trackParameters()->begin()));
1928 bool tmpgetmat = cache.m_getmaterialfromtrack;
1934 cache.m_getmaterialfromtrack =
false;
1940 trajectory.trackStates().reserve(inputTrack.trackStateOnSurfaces()->size());
1942 const Surface *firsthitsurf =
nullptr;
1943 const Surface *lasthitsurf =
nullptr;
1945 bool hasmuon =
false;
1946 bool iscombined =
false;
1947 bool seenphimeas =
false;
1951 for (; itStates != endState; ++itStates) {
1952 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1954 (pMeasurement ==
nullptr) &&
1955 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1956 matEffects != inputTrack.info().particleHypothesis()
1961 if (pMeasurement !=
nullptr) {
1962 const Surface *surf = &pMeasurement->associatedSurface();
1963 Identifier hitid = surf->associatedDetectorElementIdentifier();
1967 hitid = crot->rioOnTrack(0).identify();
1971 if (firsthitsurf ==
nullptr) {
1972 firsthitsurf = surf;
1980 if ((**itStates).trackParameters() !=
nullptr) {
1981 lastidpar = (**itStates).trackParameters();
1982 if (firstidpar ==
nullptr) {
1983 firstidpar = lastidpar;
1988 Amg::Vector3D measdir = surf->transform().rotation().col(0);
1990 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1991 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1993 if (std::abs(surf->center().z()) > 13000) {
1996 if (surf->center().perp() > 9000 && std::abs(surf->center().z()) < 13000) {
2008 if (iscombined && seenphimeas && (phiem || phibo)) {
2018 cache.m_getmaterialfromtrack && trajectory.numberOfScatterers() != 0 &&
2019 (hasmuon || cache.m_acceleration)
2021 cache.m_matfilled =
true;
2024 if (firstidpar == lastidpar) {
2025 firstidpar = lastidpar =
nullptr;
2030 !cache.m_matfilled &&
2032 m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier()) !=
2035 (firstidpar !=
nullptr)
2037 if (
m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier())) {
2044 bool tmpacc = cache.m_acceleration;
2046 bool tmpsirecal = cache.m_sirecal;
2047 std::unique_ptr<Track> tmptrack =
nullptr;
2051 cache.m_fiteloss =
true;
2052 cache.m_sirecal =
false;
2055 cache.m_asymeloss =
true;
2058 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2059 cache.m_sirecal = tmpsirecal;
2061 if (tmptrack ==
nullptr) {
2062 cache.m_matfilled =
false;
2063 cache.m_getmaterialfromtrack = tmpgetmat;
2064 cache.m_acceleration = tmpacc;
2065 cache.m_fiteloss = tmpfiteloss;
2070 bool isbrem =
false;
2072 unsigned int n_brem=0;
2074 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
2075 GXFMaterialEffects *meff = state->materialEffects();
2077 if (meff !=
nullptr) {
2081 const MaterialProperties *matprop = meff->materialProperties();
2083 double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->delta_p());
2085 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->position()));
2086 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2087 double costracksurf = 1.;
2089 costracksurf = std::abs(layerNormal.dot(layerpars->momentum().unit()));
2091 double oldde = meff->deltaE();
2093 std::unique_ptr<EnergyLoss> eloss;
2094 double sigmascat = 0;
2096 if (matprop !=
nullptr) {
2097 eloss = std::make_unique<EnergyLoss>(
2098 m_elosstool->energyLoss(*matprop,
p, 1. / costracksurf,
2100 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop,
p, 1. / costracksurf, matEffects));
2102 if (eloss !=
nullptr) {
2103 meff->setDeltaE(eloss->deltaE());
2106 MaterialProperties tmpprop(1., meff->x0(), 0., 0., 0., 0.);
2107 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop,
p, 1. / costracksurf, matEffects));
2110 meff->setScatteringSigmas(
2118 meff->setDeltaE(oldde);
2119 if (!meff->isKink()) {
2120 meff->setSigmaDeltaE(0);
2123 bremdp = meff->delta_p();
2125 }
else if (eloss !=
nullptr) {
2126 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
2128 if (meff->sigmaDeltaE() > 0) {
2134 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2135 minpar=trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2136 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2140 cache.m_matfilled =
true;
2144 trajectory.brems().clear();
2146 trajectory.brems().resize(1);
2147 trajectory.brems()[0] = bremdp;
2150 cache.m_asymeloss =
false;
2151 trajectory.setNumberOfScatterers(nscats);
2154 trajectory.setNumberOfBrems(n_brem);
2158 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2160 bool pseudoupdated =
false;
2162 if ((
track !=
nullptr) && hasid && hasmuon) {
2163 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
2165 (pseudostate ==
nullptr) ||
2167 pseudostate->fitQuality().chiSquared() < 10
2173 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2175 pseudostate->measurement()->localParameters(),
2176 pseudostate->measurement()->localCovariance()
2179 if (updpar ==
nullptr) {
2186 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2189 pseudopar->associatedSurface()
2192 pseudostate->setMeasurement(std::move(newpseudo));
2196 pseudostate->setMeasurementErrors(
errors);
2197 pseudoupdated =
true;
2200 if (pseudoupdated) {
2201 trajectory.setConverged(
false);
2202 cache.m_matfilled =
true;
2204 cache.m_matfilled =
false;
2208 cache.m_matfilled =
false;
2209 cache.m_getmaterialfromtrack = tmpgetmat;
2210 cache.m_acceleration = tmpacc;
2211 cache.m_fiteloss = tmpfiteloss;
2213 if (
track !=
nullptr) {
2215 const TrackInfo& old_info = inputTrack.info();
2216 track->info().addPatternReco(old_info);
2219 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 7477 of file GlobalChi2Fitter.cxx.
7487 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7488 ctx,
src, dst.associatedSurface(), propdir,
false
7500 &
rv.front()->associatedSurface() == &dst.associatedSurface() ||
7501 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7502 trackParametersClose(*
rv.front(),
src, 0.001) ||
7503 trackParametersClose(*
rv.front(), *dst.trackParameters(), 0.001)
7506 rv.front().reset(
nullptr);
7516 &
rv.back()->associatedSurface() == &dst.associatedSurface() ||
7517 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7518 trackParametersClose(*
rv.back(),
src, 0.001) ||
7519 trackParametersClose(*
rv.back(), *dst.trackParameters(), 0.001)
7522 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 7028 of file GlobalChi2Fitter.cxx.
7041 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7047 if (
tp ==
nullptr) {
7056 const TrkDetElementBase * de =
tp->associatedSurface().associatedDetectorElement();
7058 if (de ==
nullptr) {
7069 if (id_set.find(
id) != id_set.end()) {
7118 if (sct_set.find(
os) != sct_set.end()) {
7119 ++
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 7213 of file GlobalChi2Fitter.cxx.
7227 constexpr
uint min_meas = 3;
7228 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7232 bool seen_meas =
false;
7234 std::set<Identifier> id_set;
7235 std::set<Identifier> sct_set;
7241 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7264 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7271 if (seen_meas && dist >= 2.5) {
7278 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7279 std::vector<std::unique_ptr<TrackParameters>>
states;
7287 if (hc.has_value()) {
7308 GXFTrackState & last =
states.back();
7319 last.trackParameters() !=
nullptr &&
7326 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7328 *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 7138 of file GlobalChi2Fitter.cxx.
7146 GXFTrackState * lastmeas =
nullptr;
7148 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7160 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7161 rv.reserve(trajectory.trackStates().size());
7168 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7182 rv.emplace_back(*
s);
7189 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7191 if (de !=
nullptr) {
7204 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 8344 of file GlobalChi2Fitter.cxx.
8354 if (cond_obj ==
nullptr) {
8359 cond_obj->getInitializedCache(cache.m_field_cache);
◆ initialize()
StatusCode Trk::GlobalChi2Fitter::initialize |
( |
| ) |
|
|
overridevirtual |
Definition at line 207 of file GlobalChi2Fitter.cxx.
227 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
228 return StatusCode::FAILURE;
253 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
277 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
278 return StatusCode::FAILURE;
283 return StatusCode::SUCCESS;
◆ isMuonTrack()
bool Trk::GlobalChi2Fitter::isMuonTrack |
( |
const Track & |
intrk1 | ) |
const |
|
private |
Definition at line 8302 of file GlobalChi2Fitter.cxx.
8303 const auto *pDataVector = intrk1.measurementsOnTrack();
8304 auto nmeas1 = pDataVector->size();
8305 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8313 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8315 if (lastMeasIsCompetingRIO){
8317 testrot = &testcrot->rioOnTrack(0);
8321 if (testrot ==
nullptr) {
8322 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8325 if(penultimateIsRIO){
8326 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8328 if (penultimateIsCompetingRIO){
8330 testrot = &testcrot->rioOnTrack(0);
8337 (testrot !=
nullptr) &&
◆ iterationsOfLastFit()
int Trk::GlobalChi2Fitter::iterationsOfLastFit |
( |
| ) |
const |
|
privatevirtual |
◆ mainCombinationStrategy()
Definition at line 587 of file GlobalChi2Fitter.cxx.
595 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
600 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
601 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
603 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
604 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
605 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
607 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
613 muontrack->trackStateOnSurfaces()->end() :
614 muontrack->trackStateOnSurfaces()->begin();
620 const MeasurementBase *closestmuonmeas =
nullptr;
621 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
623 while (closestmuonmeas ==
nullptr) {
624 closestmuonmeas =
nullptr;
627 if ((**tsosit).measurementOnTrack() !=
nullptr) {
628 closestmuonmeas = (**tsosit).measurementOnTrack();
630 if (thispar !=
nullptr) {
631 const AmgVector(5) & parvec = thispar->parameters();
632 tp_closestmuon=thispar->associatedSurface().createUniqueTrackParameters(
633 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
647 std::unique_ptr<const TrackParameters> tmppar;
649 if (cache.m_msEntrance ==
nullptr) {
653 cache.m_msEntrance =
geometry->trackingVolume(
"MuonSpectrometerEntrance");
656 if (cache.m_msEntrance ==
nullptr) {
661 if ((tp_closestmuon !=
nullptr) && (cache.m_msEntrance !=
nullptr)) {
663 ctx, *tp_closestmuon, *cache.m_msEntrance, propdir,
nonInteracting);
666 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
668 if (tmppar !=
nullptr) {
669 const Surface & associatedSurface = tmppar->associatedSurface();
670 std::unique_ptr<Surface> muonsurf =
nullptr;
674 const CylinderBounds *cylbounds =
static_cast <const CylinderBounds *
>(&associatedSurface.bounds());
676 double radius = cylbounds->r();
677 double hlength = cylbounds->halflengthZ();
678 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
683 associatedSurface.center().z() > 0 ?
684 associatedSurface.center().z() + 1 :
685 associatedSurface.center().z() - 1
689 associatedSurface.center().x(),
690 associatedSurface.center().y(),
694 trans.translation() << newpos;
696 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&associatedSurface.bounds());
697 double rmin = discbounds->rMin();
698 double rmax = discbounds->rMax();
699 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
703 if (muonsurf !=
nullptr) {
715 std::vector<const TrackStateOnSurface *> tmp_matvec;
717 if ((matvec !=
nullptr) && !matvec->empty()) {
718 tmp_matvec = *matvec;
719 delete tmp_matvec.back();
720 tmp_matvec.pop_back();
722 for (
auto &
i : tmp_matvec) {
727 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(
i->materialEffectsOnTrack());
729 const Surface *matsurf = &meff->associatedSurface();
736 trajectory.m_fieldprop,
741 if (tmppar ==
nullptr) {
749 trajectory.m_fieldprop,
755 if (tmppar ==
nullptr) {
763 double de = std::abs(meff->energyLoss()->deltaE());
765 double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
772 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
773 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
789 for (; itStates != endState; ++itStates) {
794 bool tmpgetmat = cache.m_getmaterialfromtrack;
796 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
798 cache.m_extmat =
false;
800 cache.m_idmat =
false;
803 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
807 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
808 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
809 cache.m_getmaterialfromtrack =
true;
815 cache.m_getmaterialfromtrack = tmpgetmat;
822 trajectory.trackStates().back()->setTrackParameters(
nullptr);
825 std::unique_ptr<const TrackParameters> firstscatpar;
826 std::unique_ptr<const TrackParameters> lastscatpar;
829 double newqoverpid = 0;
832 double de = std::abs(calomeots[1].energyLoss()->deltaE());
833 double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
835 double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
836 double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
837 double elosspull = (pbefore - pafter - de) / sigmade;
839 if (std::abs(elosspull) > 10) {
840 if (elosspull > 10) {
841 newqoverpid = 1 / (de + pafter + 10 * sigmade);
843 newqoverpid = 1 / (de + pafter - 10 * sigmade);
846 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
850 const AmgVector(5) & newpar = firstidpar->parameters();
851 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
852 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
860 if (lastidpar ==
nullptr) {
866 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
867 calomeots[0].associatedSurface(),
870 trajectory.m_fieldprop,
874 if (firstscatpar ==
nullptr) {
880 *(firstismuon ? firstidpar : tp_closestmuon),
881 calomeots[2].associatedSurface(),
884 trajectory.m_fieldprop,
888 if (lastscatpar ==
nullptr) {
892 std::optional<TransportJacobian> jac1, jac2;
893 std::unique_ptr<const TrackParameters> elosspar;
895 double firstscatphi = 0;
896 double secondscatphi = 0;
897 double firstscattheta = 0;
898 double secondscattheta = 0;
899 double muonscatphi = 0;
900 double muonscattheta = 0;
902 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
903 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
905 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
907 Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
908 muonscatphi = calosegment.phi() - muonscatpar->parameters()[
Trk::phi];
910 if (std::abs(std::abs(muonscatphi) - 2 *
M_PI) < std::abs(muonscatphi)) {
911 muonscatphi += (muonscatphi < 0 ? 2 *
M_PI : -2 *
M_PI);
914 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
915 std::unique_ptr<const TrackParameters> startPar =
unique_clone(cache.m_idmat ? lastidpar.get() : indettrack->perigeeParameters());
917 for (
int i = 0;
i < 2;
i++) {
918 std::unique_ptr<const TrackParameters> tmpelosspar;
920 params1[
Trk::
phi] += muonscatphi;
921 params1[
Trk::
theta] += muonscattheta;
927 std::unique_ptr<const TrackParameters> tmppar1(muonscatpar->associatedSurface().createUniqueTrackParameters(
928 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
940 trajectory.m_fieldprop,
949 calomeots[1].associatedSurface(),
951 trajectory.m_fieldprop
955 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
964 const AmgVector(5) & newpars = tmpelosspar->parameters();
965 std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
966 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
970 elosspar = std::move(tmpelosspar);
973 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
977 calomeots[0].associatedSurface() :
978 calomeots[2].associatedSurface(),
981 trajectory.m_fieldprop,
991 calomeots[0].associatedSurface() :
992 calomeots[2].associatedSurface(),
996 trajectory.m_fieldprop
1000 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
1005 for (
int j = 0; j < 5; j++) {
1006 for (
int k = 0;
k < 5;
k++) {
1008 for (
int l = 0;
l < 5;
l++) {
1009 jac3[j][
k] += (*jac2) (j,
l) * (*jac1) (
l,
k);
1018 jac4(0, 0) = jac3[0][2];
1019 jac4(1, 1) = jac3[1][3];
1020 jac4(0, 1) = jac3[0][3];
1021 jac4(1, 0) = jac3[1][2];
1023 jac4 = jac4.inverse();
1035 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1037 if (cylsurf !=
nullptr) {
1039 if (std::abs(std::abs(dloc1) -
length) < std::abs(dloc1)) {
1048 if (discsurf !=
nullptr) {
1049 if (std::abs(std::abs(dloc2) - 2 *
M_PI) < std::abs(dloc2)) {
1058 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1059 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1065 muonscatphi += dphi;
1066 muonscattheta += dtheta;
1068 double idscatphi = idscatpar->parameters()[
Trk::phi] - (scat2->parameters()[
Trk::phi] + dphi);
1070 if (std::abs(std::abs(idscatphi) - 2 *
M_PI) < std::abs(idscatphi)) {
1071 idscatphi += ((idscatphi < 0) ? 2 *
M_PI : -2 *
M_PI);
1074 double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1077 firstscatphi = muonscatphi;
1078 secondscatphi = idscatphi;
1079 firstscattheta = muonscattheta;
1080 secondscattheta = idscattheta;
1082 firstscatphi = -idscatphi;
1083 secondscatphi = -muonscatphi;
1084 firstscattheta = -idscattheta;
1085 secondscattheta = -muonscattheta;
1088 if (
i == 1 && cache.m_field_cache.toroidOn() && !firstismuon) {
1090 params2[
Trk::
phi] += idscatphi;
1091 params2[
Trk::
theta] += idscattheta;
1097 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1098 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1100 idscatpar = firstscatpar.get();
1104 *cache.m_caloEntrance,
1108 if (startPar !=
nullptr) {
1117 rot.col(2) = trackdir;
1120 trans.linear().matrix() << rot;
1121 trans.translation() << startPar->position() - .1 * trackdir;
1122 PlaneSurface curvlinsurf(trans);
1132 if (curvlinpar !=
nullptr) {
1133 startPar.reset(curvlinpar);
1137 firstscatpar = std::move(scat2);
1141 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1142 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1143 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1145 double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1146 double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1149 for (
auto &
i : tmp_matvec) {
1154 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1155 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1158 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1160 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1163 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1165 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1166 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1167 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1170 for (
auto &
i : tmp_matvec) {
1177 if (startPar ==
nullptr) {
1182 (pull1 > 5 || pull2 > 5) &&
1188 bool largegap =
false;
1189 double previousz = 0;
1191 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1192 const MaterialEffectsBase *meff = (*itStates2)->materialEffectsOnTrack();
1194 const MeasurementBase *meas = (*itStates2)->measurementOnTrack();
1196 if (meff !=
nullptr) {
1198 const MaterialEffectsOnTrack *mefot{};
1200 mefot =
static_cast<const MaterialEffectsOnTrack *
>(meff);
1202 if ( mefot and mefot->energyLoss() and
1203 std::abs(mefot->energyLoss()->deltaE()) > 250 &&
1204 mefot->energyLoss()->sigmaDeltaE() < 1.e-9
1209 cache.m_extmat =
false;
1211 cache.m_idmat =
false;
1217 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1226 itStates2 == endState2 - 1 &&
1228 tpar->position().perp() > 9000 &&
1229 std::abs(tpar->position().z()) < 13000
1231 std::unique_ptr<const TrackParameters> pseudopar(tpar->clone());
1235 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1238 pseudopar->associatedSurface()
1241 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1248 pseudostate->setMeasurementErrors(
errors);
1249 trajectory.addMeasurementState(std::move(pseudostate));
1255 std::abs(trajectory.trackStates().back()->position().z()) > 20000 &&
1256 std::abs(previousz) < 12000
1262 previousz = trajectory.trackStates().back()->position().z();
1272 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
◆ makePerigee()
Definition at line 4588 of file GlobalChi2Fitter.cxx.
4593 const PerigeeSurface *persurf =
nullptr;
4596 persurf =
static_cast<const PerigeeSurface *
>(¶m.associatedSurface());
4598 if ((persurf !=
nullptr) && (!cache.m_acceleration || persurf->center().perp() > 5)) {
4600 return param.associatedSurface().createUniqueTrackParameters(
4605 if (cache.m_acceleration) {
4609 PerigeeSurface tmppersf;
4610 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4611 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4613 if (per ==
nullptr) {
4614 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
◆ makeProtoState()
Definition at line 2601 of file GlobalChi2Fitter.cxx.
2613 ) && cache.m_getmaterialfromtrack
2615 if (cache.m_acceleration && trajectory.numberOfHits() == 0) {
2621 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(tsos->materialEffectsOnTrack());
2623 std::unique_ptr<GXFMaterialEffects> newmeff;
2626 meff->scatteringAngles() or
2627 meff->energyLoss() or
2629 (tsos->trackParameters() ==
nullptr)
2631 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2635 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2637 std::abs(1. / tsos->trackParameters()->parameters()[
Trk::qOverP]),
2650 meff->thicknessInX0(),
2656 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2660 (meff->energyLoss() !=
nullptr) &&
2661 meff->energyLoss()->sigmaDeltaE() > 0 &&
2664 ((meff->scatteringAngles() ==
nullptr) || meff->thicknessInX0() == 0)
2667 newmeff->setSigmaDeltaE(meff->energyLoss()->sigmaDeltaE());
2670 (tsos->trackParameters() !=
nullptr) &&
2671 !trajectory.trackStates().empty() &&
2672 ((**trajectory.trackStates().rbegin()).trackParameters() !=
nullptr)
2674 double delta_p = 1000 * (
2675 tsos->trackParameters()->parameters()[
Trk::qOverP] -
2676 (**trajectory.trackStates().rbegin()).trackParameters()->
2680 newmeff->setdelta_p(delta_p);
2684 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(newmeff),
unique_clone(tsos->trackParameters())),
index);
2691 bool isoutlier =
false;
2700 tsos->measurementOnTrack(),
2701 tsos->trackParameters(),
◆ makeProtoStateFromMeasurement()
Definition at line 2708 of file GlobalChi2Fitter.cxx.
2718 if (!measbase->associatedSurface().associatedDetectorElementIdentifier().is_valid()) {
2720 seg =
static_cast<const Segment *
>(measbase);
2727 imax = (
int) seg->numberOfMeasurementBases();
2730 for (
int i = 0;
i <
imax;
i++) {
2731 const MeasurementBase *measbase2 = ((seg !=
nullptr) &&
m_decomposesegments) ? seg->measurement(
i) : measbase;
2733 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2734 std::unique_ptr<const MeasurementBase>(measbase2->clone()),
2735 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->clone() :
nullptr)
2738 double sinstereo = 0;
2742 Identifier hitid = measbase2->associatedSurface().associatedDetectorElementIdentifier();
2747 hitid = crot->rioOnTrack(0).identify();
2751 bool measphi =
false;
2754 bool rotated =
false;
2770 if (measbase->localParameters().parameterKey() != 1) {
2792 const double diagonalProduct =
covmat(0, 0) *
covmat(1, 1);
2793 const double element01Sq =
covmat(0, 1) *
covmat(0, 1);
2794 const double sqrtTerm = std::sqrt(
2795 (traceCov) * (traceCov) - 4. * (diagonalProduct - element01Sq)
2801 sinstereo =
std::sin(0.5 * std::asin(2 *
covmat(0, 1) / (-sqrtTerm)));
2815 const Surface *surf = &measbase2->associatedSurface();
2816 Amg::Vector3D measdir = surf->transform().rotation().col(0);
2818 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
2819 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2827 int param_index = 0;
2834 errors[1] = std::sqrt(
covmat(param_index, param_index));
2839 errors[2] = std::sqrt(
covmat(param_index, param_index));
2844 errors[3] = std::sqrt(
covmat(param_index, param_index));
2849 errors[4] = std::sqrt(
covmat(param_index, param_index));
2854 ATH_MSG_DEBUG(
"PseudoMeasurement, pos=" << measbase2->globalPosition());
2857 ATH_MSG_DEBUG(
"VertexOnTrack, pos=" << measbase2->globalPosition());
2860 ATH_MSG_DEBUG(
"Segment, pos=" << measbase2->globalPosition());
2870 ptsos->setMeasurementErrors(
errors);
2871 ptsos->setSinStereo(sinstereo);
2872 ptsos->setMeasurementType(hittype);
2873 ptsos->setMeasuresPhi(measphi);
2875 if (isoutlier && !cache.m_reintoutl) {
2880 bool ok = trajectory.addMeasurementState(std::move(ptsos),
index);
◆ makeTrack()
Definition at line 7346 of file GlobalChi2Fitter.cxx.
7353 auto trajectory = std::make_unique<Trk::TrackStates>();
7359 GXFTrajectory tmptrajectory(oldtrajectory);
7361 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7363 if (perigee_ts ==
nullptr) {
7367 tmptrajectory.addBasicState(std::move(perigee_ts), cache.m_acceleration ? 0 : tmptrajectory.numberOfUpstreamStates());
7369 trajectory->reserve(tmptrajectory.trackStates().size());
7370 for (
auto & hit : tmptrajectory.trackStates()) {
7375 hit->resetTrackCovariance();
7381 hit->materialEffects()))
7387 auto trackState = hit->trackStateOnSurface();
7388 hit->resetTrackCovariance();
7389 trajectory->emplace_back(trackState.release());
7392 auto qual = std::make_unique<FitQuality>(tmptrajectory.chi2(), tmptrajectory.nDOF());
7403 if (matEffects ==
electron && tmptrajectory.hasKink()) {
7408 if (tmptrajectory.m_straightline) {
7412 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7422 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7431 std::optional<TrackHoleCount> hole_count;
7456 if (hole_count.has_value()) {
7469 rv->setTrackSummary(std::move(
ts));
◆ makeTrackFillDerivativeMatrix()
void Trk::GlobalChi2Fitter::makeTrackFillDerivativeMatrix |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
oldtrajectory |
|
) |
| |
|
staticprivate |
Definition at line 6810 of file GlobalChi2Fitter.cxx.
6814 Amg::MatrixX & derivs = oldtrajectory.weightedResidualDerivatives();
6818 for (
auto & hit : oldtrajectory.trackStates()) {
6819 if (
const auto *pMeas{hit->measurement()};
6825 nrealmeas += hit->numberOfMeasuredParameters();
6828 cache.m_derivmat.resize(nrealmeas, oldtrajectory.numberOfFitParameters());
6829 cache.m_derivmat.setZero();
6832 int nperpars = oldtrajectory.numberOfPerigeeParameters();
6833 int nscat = oldtrajectory.numberOfScatterers();
6834 for (
auto & hit : oldtrajectory.trackStates()) {
6835 if (
const auto *pMeas{hit->measurement()};
6841 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
6842 for (
int j = 0; j < oldtrajectory.numberOfFitParameters(); j++) {
6843 cache.m_derivmat(
i, j) = derivs(measindex2, j) *
errors[measindex2];
6844 if ((j == 4 && !oldtrajectory.m_straightline) || j >= nperpars + 2 * nscat) {
6845 cache.m_derivmat(
i, j) *= 1000;
6852 measindex += hit->numberOfMeasuredParameters();
6853 }
else if (hit->materialEffects() ==
nullptr) {
6854 measindex2 += hit->numberOfMeasuredParameters();
◆ makeTrackFindPerigee()
◆ makeTrackFindPerigeeParameters()
Definition at line 6859 of file GlobalChi2Fitter.cxx.
6865 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
6866 std::tie(firstmeasstate, lastmeasstate) = oldtrajectory.findFirstLastMeasurement();
6867 std::unique_ptr<const TrackParameters> per(
nullptr);
6870 std::unique_ptr<const TrackParameters> prevpar(
6871 firstmeasstate->trackParameters() !=
nullptr ?
6872 firstmeasstate->trackParameters()->clone() :
6875 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.upstreamMaterialLayers();
6878 for (
int i = (
int)upstreamlayers.size() - 1;
i >= 0;
i--) {
6879 if (prevpar ==
nullptr) {
6886 if (
layer ==
nullptr) {
6887 layer = upstreamlayers[
i].second;
6890 DistanceSolution distsol =
layer->surfaceRepresentation().straightLineDistanceEstimate(
6891 prevpar->position(), prevpar->momentum().unit()
6895 if (distsol.numberOfSolutions() == 2) {
6900 if (distsol.first() * distsol.second() < 0 && !
first) {
6909 std::unique_ptr<const TrackParameters> layerpar(
6913 layer->surfaceRepresentation(),
6916 oldtrajectory.m_fieldprop,
6921 if (layerpar ==
nullptr) {
6925 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
6929 prevpar = std::move(layerpar);
6933 const Layer *startlayer = firstmeasstate->trackParameters()->associatedSurface().associatedLayer();
6935 if ((startlayer !=
nullptr) && (startlayer->layerMaterialProperties() !=
nullptr)) {
6936 double startfactor = startlayer->layerMaterialProperties()->alongPostFactor();
6937 const Surface & discsurf = startlayer->surfaceRepresentation();
6940 startfactor = startlayer->layerMaterialProperties()->oppositePostFactor();
6942 if (startfactor > 0.5) {
6943 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6944 firstmeasstate->trackParameters(), *startlayer,
oppositeMomentum, matEffects
6947 if (updatedpar !=
nullptr) {
6948 firstmeasstate->setTrackParameters(std::move(updatedpar));
6957 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
6959 if ((endlayer !=
nullptr) && (endlayer->layerMaterialProperties() !=
nullptr)) {
6960 double endfactor = endlayer->layerMaterialProperties()->alongPreFactor();
6961 const Surface & discsurf = endlayer->surfaceRepresentation();
6964 endfactor = endlayer->layerMaterialProperties()->oppositePreFactor();
6967 if (endfactor > 0.5) {
6968 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6969 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
6972 if (updatedpar !=
nullptr) {
6973 lastmeasstate->setTrackParameters(std::move(updatedpar));
6978 if (prevpar !=
nullptr) {
6985 oldtrajectory.m_fieldprop,
6990 if (per ==
nullptr) {
6991 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
6996 }
else if (cache.m_acceleration && (firstmeasstate->trackParameters() !=
nullptr)) {
6998 *firstmeasstate->trackParameters(),
7004 per.reset(oldtrajectory.referenceParameters()->clone());
◆ myfit()
Definition at line 4621 of file GlobalChi2Fitter.cxx.
4629 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4633 if (!trajectory.m_straightline) {
4634 if (trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
4635 trajectory.m_straightline = !cache.m_field_cache.solenoidOn();
4636 }
else if ((trajectory.prefit() == 0) && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0) {
4637 trajectory.m_straightline = !cache.m_field_cache.toroidOn();
4639 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
4644 cache.m_lastiter = 0;
4648 if (trajectory.numberOfPerigeeParameters() == -1) {
4649 cache.incrementFitStatus(
S_FITS);
4650 if (trajectory.m_straightline) {
4651 trajectory.setNumberOfPerigeeParameters(4);
4653 trajectory.setNumberOfPerigeeParameters(5);
4657 if (trajectory.nDOF() < 0) {
4662 cache.m_phiweight.clear();
4663 cache.m_firstmeasurement.clear();
4664 cache.m_lastmeasurement.clear();
4667 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4671 if (matEffects ==
Trk::electron && trajectory.m_straightline) {
4677 trajectory.setMass(
mass);
4679 ATH_MSG_DEBUG(
"start param: " << param <<
" pos: " << param.position() <<
" pt: " << param.pT());
4681 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4683 if (!cache.m_acceleration && (per ==
nullptr)) {
4693 cache.m_acceleration &&
4694 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits() &&
4697 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4700 cache.m_fastmat =
false;
4705 !cache.m_acceleration ||
4706 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() != trajectory.numberOfHits()
4708 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4711 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4715 if (cache.m_acceleration && (trajectory.referenceParameters() ==
nullptr) && (per ==
nullptr)) {
4718 if (trajectory.numberOfScatterers() >= 2) {
4719 GXFTrackState *scatstate =
nullptr;
4720 GXFTrackState *scatstate2 =
nullptr;
4723 for (std::vector<std::unique_ptr<GXFTrackState>>::
iterator it =
4724 trajectory.trackStates().begin();
4725 it != trajectory.trackStates().end();
4729 scatindex == trajectory.numberOfScatterers() / 2 ||
4730 (**it).materialEffects()->deltaE() == 0
4732 scatstate2 = (*it).get();
4737 scatstate = (*it).get();
4744 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4745 throw std::logic_error(
"Invalid scatterer");
4748 vertex = .49 * (scatstate->position() + scatstate2->position());
4750 int nstates = (
int) trajectory.trackStates().size();
4752 trajectory.trackStates()[nstates / 2 - 1]->position() +
4753 trajectory.trackStates()[nstates / 2]->position()
4757 PerigeeSurface persurf(
vertex);
4758 std::unique_ptr<const TrackParameters> nearestpar;
4759 double mindist = 99999;
4760 std::vector < GXFTrackState * >mymatvec;
4762 for (
auto &
it : trajectory.trackStates()) {
4763 if ((*it).trackParameters() ==
nullptr) {
4768 .straightLineDistanceEstimate(
4769 (*it).trackParameters()->position(),
4770 (*it).trackParameters()->momentum().unit())
4774 (cache.m_caloEntrance ==
nullptr) ||
4775 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4779 (((*it).measurement() !=
nullptr) && insideid) || (
4780 ((*it).materialEffects() !=
nullptr) &&
4782 (*it).materialEffects()->deltaE() == 0 ||
4783 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4785 (*it).materialEffects()->deltaPhi() != 0
4789 double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4790 if (dist < mindist) {
4798 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4799 mymatvec.push_back(
it.get());
4803 if (nearestpar ==
nullptr) {
4807 for (
auto & state : mymatvec) {
4809 const Surface &matsurf = state->associatedSurface();
4810 DistanceSolution distsol = matsurf.straightLineDistanceEstimate(
4811 nearestpar->position(), nearestpar->momentum().unit());
4815 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4819 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4825 trajectory.m_fieldprop,
4829 if (tmppar ==
nullptr) {
4837 trajectory.m_fieldprop,
4841 if (tmppar ==
nullptr) {
4853 if (state->materialEffects()->sigmaDeltaE() > 0) {
4854 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4857 double de = std::abs(state->materialEffects()->deltaE());
4859 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4865 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4866 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4870 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4876 trajectory.m_fieldprop,
4881 if (tmpPars !=
nullptr) {
4882 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4887 double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4888 double toteloss = std::abs(trajectory.totalEnergyLoss());
4889 double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4893 per = per->associatedSurface().createUniqueTrackParameters(
4898 if (per ==
nullptr) {
4906 PerigeeSurface persurf2(per->position());
4907 per = persurf2.createUniqueTrackParameters(
4915 }
else if (per ==
nullptr) {
4919 if ((per ==
nullptr) && (trajectory.referenceParameters() ==
nullptr)) {
4927 if (trajectory.m_straightline && (per !=
nullptr)) {
4928 if (trajectory.numberOfPerigeeParameters() == -1) {
4929 trajectory.setNumberOfPerigeeParameters(4);
4933 per = per->associatedSurface().createUniqueTrackParameters(
4936 }
else if (trajectory.numberOfPerigeeParameters() == -1) {
4937 trajectory.setNumberOfPerigeeParameters(5);
4940 if (per !=
nullptr) {
4941 trajectory.setReferenceParameters(std::move(per));
4944 int nfitpar = trajectory.numberOfFitParameters();
4945 int nperpars = trajectory.numberOfPerigeeParameters();
4946 int nscat = trajectory.numberOfScatterers();
4947 int nbrem = trajectory.numberOfBrems();
4950 Eigen::MatrixXd a_inv;
4951 a.resize(nfitpar, nfitpar);
4956 derivPool.setZero();
4958 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
4959 if (state->materialEffects() !=
nullptr) {
4962 state->setDerivatives(derivPool);
4965 bool doderiv =
true;
4967 int tmpminiter = cache.m_miniter;
4970 cache.m_lastiter =
it;
4976 cache.m_miniter = tmpminiter;
4980 if (!trajectory.converged()) {
4981 cache.m_fittercode =
4991 cache.m_miniter = tmpminiter;
4995 int nhits = trajectory.numberOfHits();
4996 int ntrthits = trajectory.numberOfTRTHits();
4997 int nsihits = trajectory.numberOfSiliconHits();
4998 double redchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
4999 double prevredchi2 = (trajectory.nDOF() > 0) ? trajectory.prevchi2() / trajectory.nDOF() : 0;
5008 (redchi2 < prevredchi2 &&
5009 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
5010 nsihits + ntrthits == nhits
5015 if (
it != 1 || nsihits != 0 || trajectory.nDOF() <= 0 || trajectory.chi2() / trajectory.nDOF() <= 3) {
5020 cache.m_miniter = tmpminiter;
5027 int ntrtprechits = trajectory.numberOfTRTPrecHits();
5028 int ntrttubehits = trajectory.numberOfTRTTubeHits();
5030 if (ntrtprechits+ntrttubehits) {
5031 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
5033 if (phf<m_minphfcut && it>=3) {
5034 if ((ntrtprechits+ntrttubehits)>=15) {
5039 <<
" | nTRTPrecHits = " << ntrtprechits
5040 <<
" | nTRTTubeHits = " << ntrttubehits
5041 <<
" | nOutliers = "
5042 << trajectory.numberOfOutliers());
5044 if (!trajectory.converged()) {
5050 cache.m_miniter = tmpminiter;
5059 cache.m_miniter = tmpminiter;
5061 if (trajectory.prefit() == 0) {
5065 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
5066 if (lltOfW.info() == Eigen::Success) {
5070 int ncols =
a.cols();
5071 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
5072 a_inv = lltOfW.solve(weightInvAMG);
5081 GXFTrajectory *finaltrajectory = &trajectory;
5083 (runOutlier || cache.m_sirecal) &&
5084 trajectory.numberOfSiliconHits() == trajectory.numberOfHits()
5092 if (finaltrajectory != &trajectory) {
5094 delete finaltrajectory;
5103 if (!cache.m_acceleration && (finaltrajectory->prefit() == 0)) {
5104 if (nperpars == 5) {
5105 for (
int i = 0;
i <
a.cols();
i++) {
5106 a_inv(4,
i) *= .001;
5107 a_inv(
i, 4) *= .001;
5111 int scatterPos = nperpars + 2 * nscat;
5112 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
5113 for (
int i = 0;
i <
a.cols();
i++) {
5114 a_inv(scatterPos,
i) *= .001;
5115 a_inv(
i, scatterPos) *= .001;
5121 int nperparams = finaltrajectory->numberOfPerigeeParameters();
5122 for (
int i = 0;
i < nperparams;
i++) {
5123 for (
int j = 0; j < nperparams; j++) {
5124 (errmat) (j,
i) = a_inv(j,
i);
5128 if (trajectory.m_straightline) {
5129 (errmat) (4, 4) = 1
e-20;
5132 const AmgVector(5) & perpars = finaltrajectory->referenceParameters()->parameters();
5133 std::unique_ptr<const TrackParameters> measper(
5134 finaltrajectory->referenceParameters()->associatedSurface().createUniqueTrackParameters(
5135 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5139 finaltrajectory->setReferenceParameters(std::move(measper));
5141 cache.m_fullcovmat = a_inv;
5145 std::unique_ptr<Track>
track =
nullptr;
5147 if (finaltrajectory->prefit() > 0) {
5148 if (finaltrajectory != &trajectory) {
5149 delete finaltrajectory;
5154 if (finaltrajectory->numberOfOutliers() <=
m_maxoutliers || !runOutlier) {
5161 double cut = (finaltrajectory->numberOfSiliconHits() ==
5162 finaltrajectory->numberOfHits())
5168 (
track !=
nullptr) && (
5169 track->fitQuality()->numberDoF() != 0 &&
5170 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5173 track.reset(
nullptr);
5177 if (
track ==
nullptr) {
5181 if (finaltrajectory != &trajectory) {
5182 delete finaltrajectory;
5185 return track.release();
◆ myfit_helper()
◆ numericalDerivatives()
Definition at line 8111 of file GlobalChi2Fitter.cxx.
8118 ParamDefsAccessor paraccessor;
8126 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8129 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8136 const Surface & previousSurface = tmpprevpar->associatedSurface();
8140 for (
int i = 0;
i < 5;
i++) {
8143 if (thisdiscsurf &&
i == 1) {
8147 vecpluseps[paraccessor.pardef[
i]] += eps[
i];
8148 vecminuseps[paraccessor.pardef[
i]] -= eps[
i];
8149 if (thiscylsurf &&
i == 0) {
8150 if (vecpluseps[0] / previousSurface.bounds().r() >
M_PI) {
8151 vecpluseps[0] -= 2 *
M_PI * previousSurface.bounds().r();
8153 if (vecminuseps[0] / previousSurface.bounds().r() < -
M_PI) {
8154 vecminuseps[0] += 2 *
M_PI * previousSurface.bounds().r();
8157 if (thisdiscsurf &&
i == 1) {
8158 if (vecpluseps[
i] >
M_PI) {
8159 vecpluseps[
i] -= 2 *
M_PI;
8161 if (vecminuseps[
i] < -
M_PI) {
8162 vecminuseps[
i] += 2 *
M_PI;
8168 std::unique_ptr<const TrackParameters> parpluseps(
8169 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8178 std::unique_ptr<const TrackParameters> parminuseps(
8179 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8189 std::unique_ptr<const TrackParameters> newparpluseps(
8200 std::unique_ptr<const TrackParameters> newparminuseps(
8215 if (newparpluseps ==
nullptr) {
8227 if (newparminuseps ==
nullptr) {
8239 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8243 for (
int j = 0; j < 5; j++) {
8244 double diff = newparpluseps->parameters()[paraccessor.pardef[j]] -
8245 newparminuseps->parameters()[paraccessor.pardef[j]];
8246 if (cylsurf && j == 0) {
8247 double length = 2 *
M_PI * surf.bounds().r();
8256 if (discsurf && j == 1) {
8257 if (std::abs(std::abs(
diff) - 2 *
M_PI) < std::abs(
diff)) {
8266 (*jac) (j,
i) =
diff / (2 * eps[
i]);
◆ processTrkVolume()
Definition at line 2890 of file GlobalChi2Fitter.cxx.
2894 if (tvol ==
nullptr) {
2901 if (confinedLayers !=
nullptr) {
2906 for (; layerIter != layerVector.end(); ++layerIter) {
2908 if (*layerIter !=
nullptr) {
2913 if ((layIndex.
value() == 0) || ((*layerIter)->layerMaterialProperties() ==
nullptr)) {
2917 const CylinderLayer *cyllay =
nullptr;
2919 cyllay =
static_cast<const CylinderLayer *
>((*layerIter));
2921 const DiscLayer *disclay =
nullptr;
2924 disclay =
static_cast<const DiscLayer *
>((*layerIter));
2926 if (disclay !=
nullptr) {
2927 if (disclay->center().z() < 0) {
2928 cache.m_negdiscs.push_back(disclay);
2930 cache.m_posdiscs.push_back(disclay);
2932 }
else if (cyllay !=
nullptr) {
2933 cache.m_barrelcylinders.push_back(cyllay);
2943 for (
size_t ib = 0 ;
ib < bsurf.size(); ++
ib) {
2944 const Layer *
layer = bsurf[
ib]->surfaceRepresentation().materialLayer();
2946 if (
layer ==
nullptr)
continue;
2950 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2954 const CylinderSurface *cylsurf =
nullptr;
2957 cylsurf =
static_cast<const CylinderSurface *
>(&
layer->surfaceRepresentation());
2959 const DiscSurface *discsurf =
nullptr;
2962 discsurf =
static_cast<const DiscSurface *
>(&
layer->surfaceRepresentation());
2964 if (discsurf !=
nullptr) {
2966 discsurf->center().z() < 0 &&
2967 std::find(cache.m_negdiscs.begin(), cache.m_negdiscs.end(),
layer) == cache.m_negdiscs.end()
2969 cache.m_negdiscs.push_back(
layer);
2971 discsurf->center().z() > 0 &&
2972 std::find(cache.m_posdiscs.begin(), cache.m_posdiscs.end(),
layer) == cache.m_posdiscs.end()
2974 cache.m_posdiscs.push_back(
layer);
2977 (cylsurf !=
nullptr) &&
2978 std::find(cache.m_barrelcylinders.begin(), cache.m_barrelcylinders.end(),
layer) == cache.m_barrelcylinders.end()
2980 cache.m_barrelcylinders.push_back(
layer);
2983 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2990 if (confinedVolumes !=
nullptr) {
2996 for (; volIter != volIterEnd; ++volIter) {
2997 if (*volIter !=
nullptr) {
◆ retrieveTrackingGeometry()
Definition at line 912 of file GlobalChi2Fitter.h.
917 if (!handle.isValid()) {
920 return handle.cptr();
◆ runIteration()
Definition at line 5656 of file GlobalChi2Fitter.cxx.
5667 int nfitpars = trajectory.numberOfFitParameters();
5668 int nperpars = trajectory.numberOfPerigeeParameters();
5669 int scatpars = 2 * trajectory.numberOfScatterers();
5670 int nupstreamstates = trajectory.numberOfUpstreamStates();
5671 int nbrem = trajectory.numberOfBrems();
5672 double oldchi2 = trajectory.chi2();
5673 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5674 int nsihits = trajectory.numberOfSiliconHits();
5675 int ntrthits = trajectory.numberOfTRTHits();
5676 int nhits = trajectory.numberOfHits();
5678 if (cache.m_phiweight.empty()) {
5679 cache.m_phiweight.assign(trajectory.trackStates().size(), 1);
5692 double newredchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
5694 ATH_MSG_DEBUG(
"old chi2: " << oldchi2 <<
"/" << trajectory.nDOF() <<
5695 "=" << oldredchi2 <<
" new chi2: " << trajectory.chi2() <<
"/" <<
5696 trajectory.nDOF() <<
"=" << newredchi2);
5698 if (trajectory.prefit() > 0 && trajectory.converged()) {
5704 std::vector < std::pair < double, double >>&scatsigmas = trajectory.scatteringSigmas();
5706 int nmeas = (
int)
res.size();
5708 const Amg::MatrixX & weight_deriv = trajectory.weightedResidualDerivatives();
5715 if (cache.m_firstmeasurement.empty()) {
5716 cache.m_firstmeasurement.resize(nfitpars);
5717 cache.m_lastmeasurement.resize(nfitpars);
5718 for (
int i = 0;
i < nperpars;
i++) {
5719 cache.m_firstmeasurement[
i] = 0;
5720 cache.m_lastmeasurement[
i] = nmeas - nbrem;
5725 for (
int i = 0;
i < (
int) trajectory.trackStates().size();
i++) {
5726 std::unique_ptr<GXFTrackState> & state = trajectory.trackStates()[
i];
5727 GXFMaterialEffects *meff = state->materialEffects();
5728 if (meff ==
nullptr) {
5729 measno += state->numberOfMeasuredParameters();
5731 if (meff !=
nullptr) {
5732 if (meff->sigmaDeltaTheta() != 0
5733 && ((trajectory.prefit() == 0) || meff->deltaE() == 0)) {
5734 int scatterPos = nperpars + 2 * scatno;
5735 if (
i < nupstreamstates) {
5736 cache.m_lastmeasurement[scatterPos] =
5737 cache.m_lastmeasurement[scatterPos + 1] = measno;
5738 cache.m_firstmeasurement[scatterPos] =
5739 cache.m_firstmeasurement[scatterPos + 1] = 0;
5741 cache.m_lastmeasurement[scatterPos] =
5742 cache.m_lastmeasurement[scatterPos + 1] = nmeas - nbrem;
5743 cache.m_firstmeasurement[scatterPos] =
5744 cache.m_firstmeasurement[scatterPos + 1] = measno;
5748 if (meff->sigmaDeltaE() > 0) {
5749 if (
i < nupstreamstates) {
5750 cache.m_firstmeasurement[nperpars + scatpars + bremno] = 0;
5751 cache.m_lastmeasurement[nperpars + scatpars + bremno] = measno;
5753 cache.m_firstmeasurement[nperpars + scatpars + bremno] = measno;
5754 cache.m_lastmeasurement[nperpars + scatpars + bremno] =
5764 if (
a.cols() != nfitpars) {
5768 for (
int k = 0;
k < nfitpars;
k++) {
5770 int maxmeas = nmeas - nbrem;
5771 maxmeas = cache.m_lastmeasurement[
k];
5772 minmeas = cache.m_firstmeasurement[
k];
5774 for (measno = minmeas; measno < maxmeas; measno++) {
5776 res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5780 if (
k == 4 ||
k >= nperpars + scatpars) {
5781 for (measno = nmeas - nbrem; measno < nmeas; measno++) {
5782 b[
k] +=
res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5787 for (
int l =
k;
l < nfitpars;
l++) {
5789 std::min(cache.m_lastmeasurement[
k], cache.m_lastmeasurement[
l]);
5792 cache.m_firstmeasurement[
l]);
5794 for (measno = minmeas; measno < maxmeas; measno++) {
5795 tmp += weight_deriv(measno,
k) * weight_deriv(measno,
l);
5797 a.fillSymmetric(
l,
k,
tmp);
5805 for (
int k = nperpars;
k < nperpars + scatpars;
k += 2) {
5806 a(
k,
k) += 1. / (scatsigmas[scatno].first * scatsigmas[scatno].first);
5807 a(
k + 1,
k + 1) += 1. / (scatsigmas[scatno].second * scatsigmas[scatno].second);
5811 for (
int measno = nmeas - nbrem; measno < nmeas; measno++) {
5812 for (
int k = 4;
k < nfitpars;
k++) {
5814 k = nperpars + scatpars;
5817 for (
int l =
k;
l < nfitpars;
l++) {
5819 l = nperpars + scatpars;
5821 double tmp =
a(
l,
k) + weight_deriv(measno,
k) * weight_deriv(measno,
l);
5822 a.fillSymmetric(
l,
k,
tmp);
5828 unsigned int scatno = 0;
5829 bool weightchanged =
false;
5831 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5832 GXFMaterialEffects *meff = thisstate->materialEffects();
5834 if (meff !=
nullptr) {
5835 const PlaneSurface *plsurf =
nullptr;
5838 plsurf =
static_cast < const PlaneSurface *
>(&thisstate->associatedSurface());
5839 if (meff->deltaE() == 0 || ((trajectory.prefit() == 0) && (plsurf !=
nullptr))) {
5840 weightchanged =
true;
5842 if (
a.cols() != nfitpars) {
5846 int scatNoIndex = 2 * scatno + nperpars;
5848 if (trajectory.prefit() == 0) {
5849 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5850 if (scatno >= cache.m_phiweight.size()) {
5852 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5853 throw std::range_error(
message.str());
5857 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5861 cache.m_phiweight[scatno] = 1.00000001;
5862 }
else if (
it == 1) {
5863 cache.m_phiweight[scatno] = 1.0000001;
5864 }
else if (
it <= 3) {
5865 cache.m_phiweight[scatno] = 1.0001;
5866 }
else if (
it <= 6) {
5867 cache.m_phiweight[scatno] = 1.01;
5869 cache.m_phiweight[scatno] = 1.1;
5872 a(scatNoIndex, scatNoIndex) *= cache.m_phiweight[scatno];
5876 else if (trajectory.prefit() >= 2) {
5877 if (newredchi2 > oldredchi2 - 1 && newredchi2 < oldredchi2) {
5878 a(scatNoIndex, scatNoIndex) *= 1.0001;
5879 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5880 }
else if (newredchi2 > oldredchi2 - 25 && newredchi2 < oldredchi2) {
5881 a(scatNoIndex, scatNoIndex) *= 1.001;
5882 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5884 a(scatNoIndex, scatNoIndex) *= 1.1;
5885 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.001;
5891 thisstate->materialEffects()->sigmaDeltaPhi() != 0 &&
5892 ((trajectory.prefit() == 0) || thisstate->materialEffects()->deltaE() == 0)
5900 (trajectory.prefit() == 2) &&
5902 trajectory.numberOfBrems() > 0 &&
5903 (newredchi2 < oldredchi2 - 25 || newredchi2 > oldredchi2)
5908 if (doderiv || weightchanged) {
5912 if (trajectory.converged()) {
5913 if ((trajectory.prefit() == 0) && nsihits + ntrthits != nhits) {
5914 unsigned int scatno = 0;
5916 if (
a.cols() != nfitpars) {
5920 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5921 if ((thisstate->materialEffects() !=
nullptr) && thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5922 if (scatno >= cache.m_phiweight.size()) {
5924 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5925 throw std::range_error(
message.str());
5928 const PlaneSurface *plsurf =
nullptr;
5931 plsurf =
static_cast<const PlaneSurface *
>(&thisstate->associatedSurface());
5933 if (thisstate->materialEffects()->deltaE() == 0 || (plsurf !=
nullptr)) {
5934 int scatNoIndex = 2 * scatno + nperpars;
5935 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5936 cache.m_phiweight[scatno] = 1;
5939 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5952 (newredchi2 < 2 || (newredchi2 < oldredchi2 && newredchi2 > oldredchi2 - .5)) &&
5953 (trajectory.prefit() == 0)
◆ 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 6316 of file GlobalChi2Fitter.cxx.
6325 bool trackok =
false;
6326 GXFTrajectory *oldtrajectory = &trajectory;
6327 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6328 GXFTrajectory *newtrajectory =
nullptr;
6329 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6335 while (!trackok && oldtrajectory->nDOF() > 0) {
6337 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->trackStates();
6340 Amg::MatrixX & weightderiv = oldtrajectory->weightedResidualDerivatives();
6341 int nfitpars = oldtrajectory->numberOfFitParameters();
6342 int nhits = oldtrajectory->numberOfHits();
6343 int nsihits = oldtrajectory->numberOfSiliconHits();
6345 if (nhits != nsihits) {
6349 double maxsipull = -1;
6351 int hitno_maxsipull = -1;
6352 int measno_maxsipull = -1;
6353 int stateno_maxsipull = 0;
6354 GXFTrackState *state_maxsipull =
nullptr;
6359 int noutl = oldtrajectory->numberOfOutliers();
6365 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6366 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6372 double *
errors = state->measurementErrors();
6374 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6375 double sinstereo = state->sinStereo();
6376 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6377 double weight1 = -1;
6379 if (hitcov(0, 0) > trackcov(0, 0)) {
6380 if (sinstereo == 0) {
6384 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6385 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6396 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6399 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6402 sipull1 =
std::max(sipull1, sipull2);
6404 if (sipull1 > maxsipull) {
6405 maxsipull = sipull1;
6406 measno_maxsipull = measno;
6407 state_maxsipull = state.get();
6408 stateno_maxsipull = stateno;
6409 hitno_maxsipull = hitno;
6420 measno += state->numberOfMeasuredParameters();
6424 double maxpull = maxsipull;
6426 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6427 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6436 oldtrajectory->chi2() / oldtrajectory->nDOF() > .25 *
m_chi2cut
6438 state_maxsipull = oldtrajectory->trackStates()[stateno_maxsipull].get();
6439 const PrepRawData *prd{};
6441 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6442 prd = rot->prepRawData();
6444 std::unique_ptr < const RIO_OnTrack > broadrot;
6445 double *olderror = state_maxsipull->measurementErrors();
6447 const TrackParameters *trackpar_maxsipull = state_maxsipull->trackParameters();
6449 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6450 std::unique_ptr<const TrackParameters> trackparForCorrect(
6451 trackpar_maxsipull->associatedSurface().createUniqueTrackParameters(
6457 state_maxsipull->hasTrackCovariance()
6459 state_maxsipull->trackCovariance())
6463 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6464 double newpull = -1;
6465 double newpull1 = -1;
6466 double newpull2 = -1;
6467 double newres1 = -1;
6468 double newres2 = -1;
6469 double newsinstereo = 0;
6473 !state_maxsipull->isRecalibrated() &&
6475 oldtrajectory->chi2() / trajectory.nDOF() > .3 *
m_chi2cut &&
6483 newerror[0] = std::sqrt(
covmat(0, 0));
6485 if (state_maxsipull->sinStereo() != 0) {
6503 newerror[0] = std::sqrt(
v0);
6506 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6508 if (cosstereo != 1.) {
6510 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6511 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6514 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6517 if (newerror[0] == 0.0) {
6518 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6522 newpull1 = std::abs(newres1 / newerror[0]);
6526 newerror[1] = std::sqrt(
covmat(1, 1));
6527 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6528 newpull2 = std::abs(newres2 / newerror[1]);
6531 newpull =
std::max(newpull1, newpull2);
6537 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6539 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6540 throw std::runtime_error(
6541 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6546 newtrajectory = oldtrajectory;
6548 if (
a.cols() != nfitpars) {
6552 double oldres1 =
res[measno_maxsipull];
6553 res[measno_maxsipull] = newres1;
6554 err[measno_maxsipull] = newerror[0];
6556 for (
int i = 0;
i < nfitpars;
i++) {
6557 if (weightderiv(measno_maxsipull,
i) == 0) {
6561 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6563 for (
int j =
i; j < nfitpars; j++) {
6567 weightderiv(measno_maxsipull,
i) *
6568 weightderiv(measno_maxsipull, j) *
6569 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6573 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6577 double oldres2 =
res[measno_maxsipull + 1];
6578 res[measno_maxsipull + 1] = newres2;
6579 err[measno_maxsipull + 1] = newerror[1];
6581 for (
int i = 0;
i < nfitpars;
i++) {
6582 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6586 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6588 for (
int j =
i; j < nfitpars; j++) {
6592 weightderiv(measno_maxsipull + 1,
i) *
6593 weightderiv(measno_maxsipull + 1, j) *
6594 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6599 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6604 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6605 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6606 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6607 olderror[1] <<
" newerror_1=" << newerror[1]
6610 state_maxsipull->setMeasurement(std::move(broadrot));
6611 state_maxsipull->setSinStereo(newsinstereo);
6612 state_maxsipull->setMeasurementErrors(newerror);
6616 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6617 (oldtrajectory->chi2() / oldtrajectory->nDOF() > .3 *
m_chi2cut || noutl > 1)
6621 (oldtrajectory->nDOF() > 1 || hittype_maxsipull ==
TrackState::SCT) &&
6626 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6627 measno_maxsipull <<
" pull=" << maxsipull
6634 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6635 newtrajectory = cleanup_newtrajectory.get();
6637 if (newa.cols() != nfitpars) {
6642 Amg::MatrixX & newweightderiv = newtrajectory->weightedResidualDerivatives();
6643 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6644 throw std::runtime_error(
6645 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6649 double oldres1 =
res[measno_maxsipull];
6650 newres[measno_maxsipull] = 0;
6652 for (
int i = 0;
i < nfitpars;
i++) {
6653 if (weightderiv(measno_maxsipull,
i) == 0) {
6657 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6659 for (
int j =
i; j < nfitpars; j++) {
6663 weightderiv(measno_maxsipull,
i) *
6664 weightderiv(measno_maxsipull, j)
6668 newweightderiv(measno_maxsipull,
i) = 0;
6672 double oldres2 =
res[measno_maxsipull + 1];
6673 newres[measno_maxsipull + 1] = 0;
6675 for (
int i = 0;
i < nfitpars;
i++) {
6676 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6680 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6682 for (
int j =
i; j < nfitpars; j++) {
6683 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6690 weightderiv(measno_maxsipull + 1,
i) *
6691 weightderiv(measno_maxsipull + 1, j)
6695 newweightderiv(measno_maxsipull + 1,
i) = 0;
6699 newtrajectory->setOutlier(stateno_maxsipull);
6705 newtrajectory->setConverged(
false);
6721 if (!newtrajectory->converged()) {
6723 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6730 if (!newtrajectory->converged()) {
6739 double oldchi2 = oldtrajectory->chi2() / oldtrajectory->nDOF();
6740 double newchi2 = (newtrajectory->nDOF() > 0) ? newtrajectory->chi2() / newtrajectory->nDOF() : 0;
6743 if (newtrajectory->nDOF() != oldtrajectory->nDOF() && maxsipull > cut2) {
6744 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6746 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6751 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6752 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6760 (void)cleanup_oldtrajectory.release();
6761 return oldtrajectory;
6763 if (oldtrajectory != newtrajectory) {
6764 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6765 oldtrajectory = newtrajectory;
6772 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6773 if (lltOfW.info() == Eigen::Success) {
6777 int ncols =
a.cols();
6778 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6779 fullcov = lltOfW.solve(weightInvAMG);
6797 oldtrajectory->nDOF() > 0 &&
6798 oldtrajectory->chi2() / oldtrajectory->nDOF() >
m_chi2cut &&
6806 (void)cleanup_oldtrajectory.release();
6807 return oldtrajectory;
◆ runTrackCleanerTRT()
Definition at line 6147 of file GlobalChi2Fitter.cxx.
6159 if (
it == 1 && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
6163 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6166 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6167 int nfitpars = trajectory.numberOfFitParameters();
6169 if (
a.cols() != nfitpars) {
6173 int nperpars = trajectory.numberOfPerigeeParameters();
6174 int nscats = trajectory.numberOfScatterers();
6177 bool outlierremoved =
false;
6178 bool hitrecalibrated =
false;
6180 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6181 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6189 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6193 trajectory.setOutlier(stateno);
6194 outlierremoved =
true;
6196 double *
errors = state->measurementErrors();
6197 double olderror =
errors[0];
6199 trajectory.updateTRTHitCount(stateno, olderror);
6201 for (
int i = 0;
i < nfitpars;
i++) {
6202 if (weightderiv(measno,
i) == 0) {
6206 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6208 for (
int j =
i; j < nfitpars; j++) {
6211 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6214 weightderiv(measno,
i) = 0;
6218 }
else if (trtrecal) {
6219 double *
errors = state->measurementErrors();
6220 double olderror =
errors[0];
6222 const auto *
const thisMeasurement{state->measurement()};
6228 if (oldrot->prepRawData() !=
nullptr) {
6229 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6231 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6233 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6234 double distance = std::abs(std::abs(trackradius) - dcradius);
6236 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6237 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters()));
6238 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6239 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters()));
6242 if (newrot !=
nullptr) {
6244 hitrecalibrated =
true;
6248 if ((measno < 0) or (measno >= (
int)
res.size())) {
6249 throw std::runtime_error(
6250 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6254 double oldres =
res[measno];
6255 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6257 state->setMeasurement(std::move(newrot));
6259 trajectory.updateTRTHitCount(stateno, olderror);
6261 for (
int i = 0;
i < nfitpars;
i++) {
6262 if (weightderiv(measno,
i) == 0) {
6266 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6268 for (
int j =
i; j < nfitpars; j++) {
6272 !cache.m_phiweight.empty() &&
6275 i < nperpars + 2 * nscats &&
6276 (
i - nperpars) % 2 == 0
6278 weight = cache.m_phiweight[(
i - nperpars) / 2];
6283 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6286 weightderiv(measno,
i) *= olderror / newerror;
6289 res[measno] = newres;
6290 err[measno] = newerror;
6299 measno += state->numberOfMeasuredParameters();
6303 if (trajectory.nDOF() < 0) {
6308 if (outlierremoved || hitrecalibrated) {
6310 trajectory.setConverged(
false);
6312 cache.m_miniter =
it + 2;
◆ setMinIterations()
void Trk::GlobalChi2Fitter::setMinIterations |
( |
int |
| ) |
|
|
privatevirtual |
Definition at line 8277 of file GlobalChi2Fitter.cxx.
8279 (
"Configure the minimum number of Iterations via jobOptions");
◆ throwFailedToGetTrackingGeomtry()
void Trk::GlobalChi2Fitter::throwFailedToGetTrackingGeomtry |
( |
| ) |
const |
|
private |
◆ trackingGeometry()
Definition at line 905 of file GlobalChi2Fitter.h.
908 if (!cache.m_trackingGeometry)
910 return cache.m_trackingGeometry;
◆ updateEnergyLoss()
Definition at line 7775 of file GlobalChi2Fitter.cxx.
7788 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7792 double newqoverp = 0;
7794 if (meff.sigmaDeltaE() <= 0) {
7799 double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.deltaE()) * std::sqrt(
mass *
mass + oldp * oldp) + meff.deltaE() * meff.deltaE();
7806 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7812 return surf.createUniqueTrackParameters(
7813 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
◆ updateFitParameters()
Definition at line 5961 of file GlobalChi2Fitter.cxx.
5969 double d0 = refpar->parameters()[
Trk::d0];
5970 double z0 = refpar->parameters()[
Trk::z0];
5973 double qoverp = refpar->parameters()[
Trk::qOverP];
5974 int nscat = trajectory.numberOfScatterers();
5975 int nbrem = trajectory.numberOfBrems();
5976 int nperparams = trajectory.numberOfPerigeeParameters();
5978 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
5981 if (llt.info() == Eigen::Success) {
5987 if (trajectory.numberOfPerigeeParameters() > 0) {
5992 qoverp = (trajectory.m_straightline) ? 0 : .001 *
result[4] + qoverp;
6001 std::vector < std::pair < double, double >>&scatangles = trajectory.scatteringAngles();
6002 std::vector < double >&delta_ps = trajectory.brems();
6004 for (
int i = 0;
i < nscat;
i++) {
6005 scatangles[
i].first +=
result[2 *
i + nperparams];
6006 scatangles[
i].second +=
result[2 *
i + nperparams + 1];
6009 for (
int i = 0;
i < nbrem;
i++) {
6010 delta_ps[
i] +=
result[nperparams + 2 * nscat +
i];
6013 std::unique_ptr<const TrackParameters> newper(
6014 trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
6019 trajectory.setReferenceParameters(std::move(newper));
6020 trajectory.setScatteringAngles(scatangles);
6021 trajectory.setBrems(delta_ps);
◆ updatePixelROTs()
Update the Pixel ROT using the current trajectory/local track parameters.
Definition at line 6026 of file GlobalChi2Fitter.cxx.
6031 if ( trajectory.numberOfSiliconHits() == 0) {
6039 const EventContext &evtctx = Gaudi::Hive::currentContext();
6042 if (!splitProbContainer.isValid()) {
6046 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6049 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6050 int nfitpars = trajectory.numberOfFitParameters();
6053 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6058 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6061 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6071 const PrepRawData *prd{};
6073 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6074 prd = rot->prepRawData();
6084 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6085 if (!splitProb.isSplit()) {
6086 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6090 std::unique_ptr < const RIO_OnTrack > newrot;
6091 double *olderror = state->measurementErrors();
6094 double newerror[5] = {-1,-1,-1,-1,-1};
6095 double newres[2] = {-1,-1};
6104 newerror[0] = std::sqrt(
covmat(0, 0));
6105 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6106 newerror[1] = std::sqrt(
covmat(1, 1));
6107 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6109 if (
a.cols() != nfitpars) {
6114 for(
int k =0;
k<2;
k++ ){
6115 double oldres =
res[measno+
k];
6116 res[measno+
k] = newres[
k];
6117 err[measno+
k] = newerror[
k];
6119 for (
int i = 0;
i < nfitpars;
i++) {
6120 if (weightderiv(measno+
k,
i) == 0) {
6124 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6126 for (
int j =
i; j < nfitpars; j++) {
6130 weightderiv(measno+
k,
i) *
6131 weightderiv(measno+
k, j) *
6132 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6136 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6140 state->setMeasurement(std::move(newrot));
6141 state->setMeasurementErrors(newerror);
◆ 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 930 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 923 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)
JetConstituentVector::iterator iterator
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.
@ DeadElement
outside the element
Gaudi::Property< double > m_outlcut
@ z
global position (cartesian)
Gaudi::Property< bool > m_numderiv
FitterStatusCode runIteration(const EventContext &ctx, Cache &, GXFTrajectory &, int, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool &) const
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.
@ numberOfSCTDeadSensors
number of TRT hits
std::vector< SharedObject< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
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
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...
ToolHandle< IRIO_OnTrackCreator > m_ROTcreator
ToolHandle< IExtrapolator > m_extrapolator
ParametersT< 5, Charged, PerigeeSurface > Perigee
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.
ToolHandle< IPropagator > m_propagator
@ loc2
generic first and second local coordinate
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.
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.
ParametersBase< 5, Charged > TrackParameters
Gaudi::Property< int > m_maxoutliers
@ 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
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 ...
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)
void fillDerivatives(GXFTrajectory &traj, bool onlybrem=false) const
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.
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
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 runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int) const
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
virtual BinnedArraySpan< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
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
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
ParametersT< 5, Charged, PlaneSurface > AtaPlane
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
FitterStatusCode updateFitParameters(GXFTrajectory &, Amg::VectorX &, const Amg::SymMatrixX &) const
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the closest approach of two lines.
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.
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...
void updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &) const
Update the Pixel ROT using the current trajectory/local track parameters.
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
Gaudi::Property< bool > m_trtrecal
bool consistentSurfaces(U)
virtual MaterialEffectsDerivedType derivedType() const =0
Returns the concrete derived type.
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
std::span< T > BinnedArraySpan
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 &)
bool is_mdt(Identifier id) const
Gaudi::Property< bool > m_useCaloTG
virtual double r() const override final
This method returns the radius.
@ OutlierLogicFailure
outlier logic failed
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_clusterSplitProbContainer
void fillResiduals(const EventContext &ctx, Cache &, GXFTrajectory &, int, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool &) const
ParametersT< 5, Charged, StraightLineSurface > AtaStraightLine
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
@ Unknown
Track fitter not defined.