|
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 EventContext &evtctx) const |
| Update the Pixel ROT using the current trajectory/local track parameters. More...
|
|
GXFTrajectory * | runTrackCleanerSilicon (const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const |
|
void | runTrackCleanerTRT (Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const |
|
PropagationResult | calculateTrackParametersPropagateHelper (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method. More...
|
|
PropagationResult | calculateTrackParametersPropagate (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possible holes. More...
|
|
std::vector< std::reference_wrapper< GXFTrackState > > | holeSearchStates (GXFTrajectory &trajectory) const |
| Extracts a collection of track states which are important for hole search. More...
|
|
std::optional< GlobalChi2Fitter::TrackHoleCount > | holeSearchProcess (const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState >> &states) const |
| Conduct a hole search between a list of states, possibly reusing existing information. More...
|
|
void | holeSearchHelper (const std::vector< std::unique_ptr< TrackParameters >> &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const |
| Helper method for the hole search that does the actual counting of holes and dead modules. More...
|
|
FitterStatusCode | calculateTrackParameters (const EventContext &ctx, GXFTrajectory &, bool) const |
|
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > | updateEnergyLoss (const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const |
|
void | calculateTrackErrors (GXFTrajectory &, Amg::SymMatrixX &, bool) const |
|
std::optional< TransportJacobian > | numericalDerivatives (const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const |
|
virtual int | iterationsOfLastFit () const |
|
virtual void | setMinIterations (int) |
|
bool | isMuonTrack (const Track &) const |
|
void | initFieldCache (const EventContext &ctx, Cache &cache) const |
| Initialize a field cache inside a fit cache object. More...
|
|
void | throwFailedToGetTrackingGeomtry () const |
|
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 7829 of file GlobalChi2Fitter.cxx.
7830 int nstatesupstream = trajectory.numberOfUpstreamStates();
7831 int nscatupstream = trajectory.numberOfUpstreamScatterers();
7832 int nbremupstream = trajectory.numberOfUpstreamBrems();
7833 int nscats = trajectory.numberOfScatterers();
7834 int nperpars = trajectory.numberOfPerigeeParameters();
7835 int nfitpars = trajectory.numberOfFitParameters();
7837 using Matrix55 = Eigen::Matrix<double, 5, 5>;
7839 Matrix55 initialjac;
7840 initialjac.setZero();
7841 initialjac(4, 4) = 1;
7843 Matrix55 jacvertex(initialjac);
7845 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.numberOfScatterers(), initialjac);
7846 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.numberOfBrems(), initialjac);
7848 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7849 GXFTrackState *prevstate =
nullptr, *state =
nullptr;
7851 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
7853 for (
bool forward : {
false,
true}) {
7855 hit_begin = nstatesupstream;
7857 scatno = nscatupstream;
7858 bremno = nbremupstream;
7860 hit_begin = nstatesupstream - 1;
7862 scatno = trajectory.numberOfUpstreamScatterers() - 1;
7863 bremno = trajectory.numberOfUpstreamBrems() - 1;
7867 int hitno = hit_begin;
7868 forward ? (hitno < hit_end) : (hitno >= hit_end);
7869 hitno += (forward ? 1 : -1)
7872 state =
states[hitno].get();
7876 if (fillderivmat && state->derivatives().cols() != nfitpars) {
7877 state->derivatives().resize(5, nfitpars);
7878 state->derivatives().setZero();
7881 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
7883 if (hitno == (forward ? hit_end - 1 : 0)) {
7884 if (!fillderivmat) {
7892 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
7894 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
7895 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
7896 jacvertex(4, 4) = jac(4, 4);
7898 int jmin = 0, jmax = 0, jcnt = 0;
7899 int lp_bgn = 0, lp_end = 0;
7903 jcnt = jmax - jmin + 1;
7905 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
7908 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7910 i == scatno + (forward ? -1 : 1) &&
7911 prevstate !=
nullptr &&
7913 (!trajectory.prefit() || prevstate->materialEffects()->deltaE() == 0)
7915 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7916 jacscat[
i](4, 4) = jac(4, 4);
7918 calculateJac(jac, jacscat[
i], jmin, jmax);
7922 Eigen::MatrixXd & derivmat = state->derivatives();
7923 int scatterPos = nperpars + 2 *
i;
7925 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
7931 jcnt = jmax - jmin + 1;
7933 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
7936 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7938 i == bremno + (forward ? -1 : 1) &&
7940 prevstate->materialEffects() &&
7941 prevstate->materialEffects()->sigmaDeltaE() > 0
7943 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7944 jacbrem[
i](4, 4) = jac(4, 4);
7946 calculateJac(jac, jacbrem[
i], jmin, jmax);
7950 Eigen::MatrixXd & derivmat = state->derivatives();
7951 int scatterPos = nperpars + 2 * nscats +
i;
7953 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
7957 calculateJac(jac, jacvertex, 0, 4);
7961 Eigen::MatrixXd & derivmat = state->derivatives();
7962 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
7964 if (nperpars == 5) {
7965 derivmat.col(4).segment(0, 4) *= .001;
7966 derivmat(4, 4) = .001 * jacvertex(4, 4);
7972 (!trajectory.prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
7974 scatno += (forward ? 1 : -1);
7978 states[hitno]->materialEffects() &&
7979 states[hitno]->materialEffects()->sigmaDeltaE() > 0
7981 bremno += (forward ? 1 : -1);
7984 prevstate =
states[hitno].get();
◆ calculateTrackErrors()
Definition at line 7991 of file GlobalChi2Fitter.cxx.
7999 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
8000 int nstatesupstream = trajectory.numberOfUpstreamStates();
8002 GXFTrackState *prevstate =
nullptr;
8003 int i = nstatesupstream;
8004 for (
int j = 0; j < (
int)
states.size(); j++) {
8005 if (j < nstatesupstream) {
8012 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8013 if (stateno == 0 || stateno == nstatesupstream) {
8014 prevstate =
nullptr;
8017 std::unique_ptr<GXFTrackState> & state =
states[
index];
8018 if (state->materialEffects() !=
nullptr) {
8019 prevstate = state.get();
8023 if (!state->hasTrackCovariance()) {
8024 state->zeroTrackCovariance();
8026 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8028 if ((prevstate !=
nullptr) &&
8032 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8035 trackerrmat = jac * prevcov * jac.transpose();
8039 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8043 const MeasurementBase *measurement = state->measurement();
8044 const Amg::MatrixX & meascov = measurement->localCovariance();
8046 ParamDefsAccessor paraccessor;
8050 bool errorok =
true;
8051 for (
int i = 0;
i < 5;
i++) {
8052 if (measurement->localParameters().contains(paraccessor.pardef[
i])) {
8054 && trackerrmat(
i,
i) > meascov(j, j)) {
8056 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8057 trackerrmat(
i,
i) = meascov(j, j);
8058 for (
int k = 0;
k < 5;
k++) {
8068 for (
int i = 0;
i < 5;
i++) {
8072 for (
int j = 0; j < 5; j++) {
8079 if (trajectory.m_straightline) {
8080 trackerrmat(4, 4) = 1
e-20;
8084 state->trackParameters();
8086 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8088 if (state->hasTrackCovariance()) {
8089 trkerrmat = (state->trackCovariance());
8091 trkerrmat = std::nullopt;
8094 const AmgVector(5) & tpars = tmptrackpar->parameters();
8095 std::unique_ptr<const TrackParameters> trackpar(
8096 tmptrackpar->associatedSurface().createUniqueTrackParameters(tpars[0],
8101 std::move(trkerrmat))
8103 state->setTrackParameters(std::move(trackpar));
8104 FitQualityOnSurface fitQual{};
8106 if (errorok && trajectory.nDOF() > 0) {
8107 fitQual =
m_updator->fullStateFitQuality(
8108 *state->trackParameters(),
8109 measurement->localParameters(),
8110 measurement->localCovariance()
8113 fitQual = FitQualityOnSurface(0, state->numberOfMeasuredParameters());
8116 state->setFitQuality(fitQual);
8118 prevstate = state.get();
◆ calculateTrackParameters()
Definition at line 7605 of file GlobalChi2Fitter.cxx.
7613 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7614 int nstatesupstream = trajectory.numberOfUpstreamStates();
7615 const TrackParameters *prevtrackpar = trajectory.referenceParameters();
7616 std::unique_ptr<const TrackParameters> tmptrackpar;
7618 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7622 DistanceSolution distsol = surf1.straightLineDistanceEstimate(
7623 prevtrackpar->position(), prevtrackpar->momentum().unit()
7630 distsol.numberOfSolutions() > 0 &&
7631 prevtrackpar != trajectory.referenceParameters()
7641 trajectory.m_fieldprop,
7648 (
rv.m_parameters !=
nullptr) &&
7649 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7655 if (
rv.m_parameters ==
nullptr) {
7656 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7657 " pos: " << prevtrackpar->position() <<
" destination surface: " << surf1);
7661 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7665 if (
rv.m_jacobian != std::nullopt) {
7667 states[hitno]->materialEffects() !=
nullptr &&
7668 states[hitno]->materialEffects()->deltaE() != 0 &&
7669 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7670 !trajectory.m_straightline
7672 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7673 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7674 double mass = trajectory.mass();
7675 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7676 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7679 states[hitno]->setJacobian(*
rv.m_jacobian);
7680 }
else if (calcderiv) {
7685 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7687 if (meff !=
nullptr && hitno != 0) {
7688 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7689 surf, *meff, *
states[hitno]->trackParameters(), trajectory.mass(), -1
7692 if (std::holds_alternative<FitterStatusCode>(
r)) {
7693 return std::get<FitterStatusCode>(
r);
7696 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7697 prevtrackpar = tmptrackpar.get();
7699 prevtrackpar = currenttrackpar;
7703 prevtrackpar = trajectory.referenceParameters();
7705 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7708 DistanceSolution distsol = surf.straightLineDistanceEstimate(prevtrackpar->position(), prevtrackpar->momentum().unit());
7712 if (
distance < 0 && distsol.numberOfSolutions() > 0 && prevtrackpar != trajectory.referenceParameters()) {
7721 trajectory.m_fieldprop,
7727 (
rv.m_parameters !=
nullptr) &&
7729 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7734 if (
rv.m_parameters ==
nullptr) {
7735 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7736 " pos: " << prevtrackpar->
7737 position() <<
" destination surface: " << surf);
7741 if (
rv.m_jacobian != std::nullopt) {
7743 states[hitno]->materialEffects() !=
nullptr &&
7744 states[hitno]->materialEffects()->deltaE() != 0 &&
7745 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7746 !trajectory.m_straightline
7748 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7749 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7750 double mass = trajectory.mass();
7751 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7754 newp = std::sqrt(newp);
7757 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7760 states[hitno]->setJacobian(*
rv.m_jacobian);
7761 }
else if (calcderiv) {
7766 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7768 if (meff !=
nullptr) {
7769 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7770 surf, *meff, *
rv.m_parameters, trajectory.mass(), +1
7773 if (std::holds_alternative<FitterStatusCode>(
r)) {
7774 return std::get<FitterStatusCode>(
r);
7777 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7780 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7781 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 7579 of file GlobalChi2Fitter.cxx.
7588 PropagationResult
rv;
7591 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7594 if (
rv.m_parameters ==
nullptr) {
7595 propdir = invertPropdir(propdir);
7598 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 7540 of file GlobalChi2Fitter.cxx.
7549 std::unique_ptr<const TrackParameters>
rv;
7550 std::optional<TransportJacobian> jac{};
7561 if (
rv !=
nullptr && calcderiv) {
7566 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7572 return PropagationResult {
7575 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 5502 of file GlobalChi2Fitter.cxx.
5508 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5512 int nscatupstream = trajectory.numberOfUpstreamScatterers();
5513 int nbremupstream = trajectory.numberOfUpstreamBrems();
5514 int nscat = trajectory.numberOfScatterers();
5515 int nbrem = trajectory.numberOfBrems();
5516 int nperparams = trajectory.numberOfPerigeeParameters();
5518 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5521 int nmeas = (
int) weightderiv.rows();
5523 ParamDefsAccessor paraccessor;
5525 for (std::unique_ptr<GXFTrackState> & state :
states) {
5528 ((state->materialEffects() ==
nullptr) || state->materialEffects()->sigmaDeltaE() <= 0)
5534 const MeasurementBase *measbase = state->measurement();
5535 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5536 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5539 double sinstereo = 0;
5542 sinstereo = state->sinStereo();
5545 double cosstereo = (sinstereo == 0) ? 1. : std::sqrt(1 - sinstereo * sinstereo);
5547 for (
int i = 0;
i < 5;
i++) {
5549 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5555 if (trajectory.numberOfPerigeeParameters() > 0) {
5556 int cols = trajectory.m_straightline ? 4 : 5;
5559 weightderiv.row(measno).head(
cols) =
5560 (derivatives.row(0).head(
cols) * cosstereo +
5561 sinstereo * derivatives.row(1).head(
cols)) /
5564 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5568 for (
int j = scatmin; j < scatmax; j++) {
5569 int index = nperparams + ((trajectory.prefit() != 1) ? 2 * j : j);
5570 double thisderiv = 0;
5574 if (
i == 0 && sinstereo != 0) {
5575 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5580 weightderiv(measno,
index) = thisderiv /
error[measno];
5582 if (trajectory.prefit() != 1) {
5585 if (
i == 0 && sinstereo != 0) {
5586 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5591 weightderiv(measno,
index) = thisderiv /
error[measno];
5595 for (
int j = bremmin; j < bremmax; j++) {
5596 double thisderiv = 0;
5597 int index = j + nperparams + 2 * nscat;
5598 if (
i == 0 && sinstereo != 0) {
5599 thisderiv = derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index);
5601 thisderiv = derivatives(
i,
index);
5603 weightderiv(measno,
index) = thisderiv /
error[measno];
5609 double *
errors = state->measurementErrors();
5610 for (
int i = 0;
i < 5;
i++) {
5617 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5622 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5624 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5625 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5627 double mass = .001 * trajectory.mass();
5629 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5633 auto multiplier = [] (
double mass,
double qOverP){
5636 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5637 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5639 if (trajectory.numberOfPerigeeParameters() > 0) {
5640 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5643 const auto bremNoBase = nperparams + 2 * nscat;
5644 if (bremno < nbremupstream) {
5645 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5646 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5647 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5650 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5651 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5652 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
◆ fillResiduals()
Definition at line 5192 of file GlobalChi2Fitter.cxx.
5204 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5209 int nbrem = trajectory.numberOfBrems();
5210 int nperpars = trajectory.numberOfPerigeeParameters();
5211 int nfitpars = trajectory.numberOfFitParameters();
5214 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5215 int nidhits = trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits();
5216 int nsihits = trajectory.numberOfSiliconHits();
5217 int ntrthits = trajectory.numberOfTRTHits();
5218 int nhits = trajectory.numberOfHits();
5219 int nmeas = (
int)
res.size();
5221 ParamDefsAccessor paraccessor;
5222 bool scatwasupdated =
false;
5224 GXFTrackState *state_maxbrempull =
nullptr;
5225 int bremno_maxbrempull = 0;
5226 double maxbrempull = 0;
5228 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5229 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5232 const MeasurementBase *measbase = state->measurement();
5239 trajectory.nDOF() != 0 &&
5240 std::abs((trajectory.prevchi2() - trajectory.chi2()) / trajectory.nDOF()) < 15 &&
5241 !state->associatedSurface().isFree() &&
5242 nidhits < trajectory.numberOfHits() &&
5243 (nperpars == 0 || nidhits > 0) &&
5244 (!state->isRecalibrated())
5249 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5252 currenttrackpar->associatedSurface()
5255 state->setMeasurement(std::move(newpseudo));
5256 measbase = state->measurement();
5259 double *
errors = state->measurementErrors();
5263 for (
int i = 0;
i < 5;
i++) {
5265 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5276 res[measno] = residuals[
i];
5278 if (
i == 2 && std::abs(std::abs(
res[measno]) - 2 *
M_PI) < std::abs(
res[measno])) {
5279 if (
res[measno] < 0) {
5288 double *
errors = state->measurementErrors();
5289 for (
int i = 0;
i < 5;
i++) {
5299 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5301 double deltaphi = state->materialEffects()->deltaPhi();
5302 double measdeltaphi = state->materialEffects()->measuredDeltaPhi();
5303 double sigmadeltaphi = state->materialEffects()->sigmaDeltaPhi();
5304 double deltatheta = state->materialEffects()->deltaTheta();
5305 double sigmadeltatheta = state->materialEffects()->sigmaDeltaTheta();
5307 if (trajectory.prefit() != 1) {
5308 b[nperpars + 2 * scatno] -= (deltaphi - measdeltaphi) / (sigmadeltaphi * sigmadeltaphi);
5309 b[nperpars + 2 * scatno + 1] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5311 b[nperpars + scatno] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5315 deltaphi * deltaphi / (sigmadeltaphi * sigmadeltaphi) +
5316 deltatheta * deltatheta / (sigmadeltatheta * sigmadeltatheta)
5322 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5323 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5325 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5326 const double pbrem = 1. / std::abs(qoverpbrem);
5327 const double p = 1. / std::abs(qoverp);
5328 const double mass = .001 * trajectory.mass();
5330 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5332 res[nmeas - nbrem + bremno] = .001 * averagenergyloss -
energy + bremEnergy;
5334 double sigde = state->materialEffects()->sigmaDeltaE();
5335 double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5336 double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5338 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5340 if (state->materialEffects()->isKink()) {
5341 maxbrempull = -999999999;
5342 state_maxbrempull =
nullptr;
5346 cache.m_asymeloss &&
5348 (trajectory.prefit() == 0) &&
5350 sigde != sigdepos &&
5353 double elosspull =
res[nmeas - nbrem + bremno] / (.001 * sigde);
5355 if (trajectory.mass() > 100) {
5356 if (elosspull < -1) {
5357 state->materialEffects()->setSigmaDeltaE(sigdepos);
5358 }
else if (elosspull > 1) {
5359 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5362 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5363 }
else if ((trajectory.numberOfTRTHits() == 0) ||
it >= 3) {
5365 !state->materialEffects()->isKink() && (
5366 (elosspull < -.2 &&
m_fixbrem == -1 && elosspull < maxbrempull) ||
5370 bremno_maxbrempull = bremno;
5371 state_maxbrempull = state.get();
5372 maxbrempull = elosspull;
5381 (trajectory.prefit() == 0) &&
5382 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5383 state->materialEffects()->isMeasuredEloss() &&
5384 res[nmeas - nbrem + bremno] / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5387 trajectory.prefit() != 0 ? trajectory.referenceParameters()
5388 :
states[hitno - 2]->trackParameters();
5391 std::vector<MaterialEffectsOnTrack> calomeots =
5396 parforcalo->associatedSurface(),
5400 if (calomeots.size() == 3) {
5401 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5402 double newres = .001 * averagenergyloss -
energy + bremEnergy;
5403 double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5405 if (std::abs(newres / newerr) < std::abs(
res[nmeas - nbrem + bremno] /
error[nmeas - nbrem + bremno])) {
5406 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5408 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5409 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5410 res[nmeas - nbrem + bremno] = newres;
5411 error[nmeas - nbrem + bremno] = newerr;
5415 state->materialEffects()->setMeasuredEloss(
false);
5423 for (; measno < nmeas; measno++) {
5424 if (
error[measno] == 0) {
5431 if (!doderiv && (scatwasupdated)) {
5435 double oldchi2 = trajectory.chi2();
5436 trajectory.setPrevChi2(oldchi2);
5437 trajectory.setChi2(
chi2);
5439 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5440 double newredchi2 = (trajectory.nDOF() > 0) ?
chi2 / trajectory.nDOF() : 0;
5443 trajectory.prefit() > 0 && (
5444 (newredchi2 < 2 &&
it != 0) ||
5445 (newredchi2 < oldredchi2 + .1 && std::abs(newredchi2 - oldredchi2) < 1 &&
it != 1)
5448 trajectory.setConverged(
true);
5451 double maxdiff = (nsihits != 0 && nsihits + ntrthits == nhits &&
chi2 < oldchi2) ? 200 : 1.;
5453 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5455 if (miniter < cache.m_miniter) {
5456 miniter = cache.m_miniter;
5460 trajectory.setConverged(
true);
5463 if ((state_maxbrempull !=
nullptr) && trajectory.converged()) {
5464 state_maxbrempull->materialEffects()->setSigmaDeltaE(
5465 10 * state_maxbrempull->materialEffects()->sigmaDeltaEPos()
5468 state_maxbrempull->materialEffects()->setKink(
true);
5469 trajectory.setConverged(
false);
5471 double olderror =
error[nmeas - nbrem + bremno_maxbrempull];
5472 double newerror = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5473 error[nmeas - nbrem + bremno_maxbrempull] = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5475 if (
a.cols() != nfitpars) {
5479 for (
int i = 0;
i < nfitpars;
i++) {
5480 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5484 for (
int j =
i; j < nfitpars; j++) {
5488 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull, j) *
5490 (1 - olderror * olderror / (newerror * newerror))
5494 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= olderror / newerror;
5497 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, ctx);
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 7489 of file GlobalChi2Fitter.cxx.
7499 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7500 ctx,
src, dst.associatedSurface(), propdir,
false
7512 &
rv.front()->associatedSurface() == &dst.associatedSurface() ||
7513 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7514 trackParametersClose(*
rv.front(),
src, 0.001) ||
7515 trackParametersClose(*
rv.front(), *dst.trackParameters(), 0.001)
7518 rv.front().reset(
nullptr);
7528 &
rv.back()->associatedSurface() == &dst.associatedSurface() ||
7529 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7530 trackParametersClose(*
rv.back(),
src, 0.001) ||
7531 trackParametersClose(*
rv.back(), *dst.trackParameters(), 0.001)
7534 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 7032 of file GlobalChi2Fitter.cxx.
7045 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7051 if (
tp ==
nullptr) {
7060 const TrkDetElementBase * de =
tp->associatedSurface().associatedDetectorElement();
7062 if (de ==
nullptr) {
7073 if (id_set.find(
id) != id_set.end()) {
7122 if (sct_set.find(
os) != sct_set.end()) {
7123 ++
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 7217 of file GlobalChi2Fitter.cxx.
7231 constexpr
uint min_meas = 3;
7232 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7236 bool seen_meas =
false;
7238 std::set<Identifier> id_set;
7239 std::set<Identifier> sct_set;
7245 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7268 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7270 bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7272 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7273 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7283 if (seen_meas && dist >= 2.5 && zStartValid) {
7290 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7291 std::vector<std::unique_ptr<TrackParameters>>
states;
7299 if (hc.has_value()) {
7320 GXFTrackState & last =
states.back();
7331 last.trackParameters() !=
nullptr &&
7338 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7340 *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 7142 of file GlobalChi2Fitter.cxx.
7150 GXFTrackState * lastmeas =
nullptr;
7152 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7164 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7165 rv.reserve(trajectory.trackStates().size());
7172 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7186 rv.emplace_back(*
s);
7193 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7195 if (de !=
nullptr) {
7208 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 8356 of file GlobalChi2Fitter.cxx.
8366 if (cond_obj ==
nullptr) {
8371 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 8314 of file GlobalChi2Fitter.cxx.
8315 const auto *pDataVector = intrk1.measurementsOnTrack();
8316 auto nmeas1 = pDataVector->size();
8317 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8325 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8327 if (lastMeasIsCompetingRIO){
8329 testrot = &testcrot->rioOnTrack(0);
8333 if (testrot ==
nullptr) {
8334 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8337 if(penultimateIsRIO){
8338 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8340 if (penultimateIsCompetingRIO){
8342 testrot = &testcrot->rioOnTrack(0);
8349 (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");
4618 if(std::abs(per->position().z())>5000.) {
4619 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
◆ 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 7358 of file GlobalChi2Fitter.cxx.
7365 auto trajectory = std::make_unique<Trk::TrackStates>();
7371 GXFTrajectory tmptrajectory(oldtrajectory);
7373 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7375 if (perigee_ts ==
nullptr) {
7379 tmptrajectory.addBasicState(std::move(perigee_ts), cache.m_acceleration ? 0 : tmptrajectory.numberOfUpstreamStates());
7381 trajectory->reserve(tmptrajectory.trackStates().size());
7382 for (
auto & hit : tmptrajectory.trackStates()) {
7387 hit->resetTrackCovariance();
7393 hit->materialEffects()))
7399 auto trackState = hit->trackStateOnSurface();
7400 hit->resetTrackCovariance();
7401 trajectory->emplace_back(trackState.release());
7404 auto qual = std::make_unique<FitQuality>(tmptrajectory.chi2(), tmptrajectory.nDOF());
7415 if (matEffects ==
electron && tmptrajectory.hasKink()) {
7420 if (tmptrajectory.m_straightline) {
7424 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7434 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7443 std::optional<TrackHoleCount> hole_count;
7468 if (hole_count.has_value()) {
7481 rv->setTrackSummary(std::move(
ts));
◆ makeTrackFillDerivativeMatrix()
void Trk::GlobalChi2Fitter::makeTrackFillDerivativeMatrix |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
oldtrajectory |
|
) |
| |
|
staticprivate |
Definition at line 6814 of file GlobalChi2Fitter.cxx.
6818 Amg::MatrixX & derivs = oldtrajectory.weightedResidualDerivatives();
6822 for (
auto & hit : oldtrajectory.trackStates()) {
6823 if (
const auto *pMeas{hit->measurement()};
6829 nrealmeas += hit->numberOfMeasuredParameters();
6832 cache.m_derivmat.resize(nrealmeas, oldtrajectory.numberOfFitParameters());
6833 cache.m_derivmat.setZero();
6836 int nperpars = oldtrajectory.numberOfPerigeeParameters();
6837 int nscat = oldtrajectory.numberOfScatterers();
6838 for (
auto & hit : oldtrajectory.trackStates()) {
6839 if (
const auto *pMeas{hit->measurement()};
6845 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
6846 for (
int j = 0; j < oldtrajectory.numberOfFitParameters(); j++) {
6847 cache.m_derivmat(
i, j) = derivs(measindex2, j) *
errors[measindex2];
6848 if ((j == 4 && !oldtrajectory.m_straightline) || j >= nperpars + 2 * nscat) {
6849 cache.m_derivmat(
i, j) *= 1000;
6856 measindex += hit->numberOfMeasuredParameters();
6857 }
else if (hit->materialEffects() ==
nullptr) {
6858 measindex2 += hit->numberOfMeasuredParameters();
◆ makeTrackFindPerigee()
◆ makeTrackFindPerigeeParameters()
Definition at line 6863 of file GlobalChi2Fitter.cxx.
6869 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
6870 std::tie(firstmeasstate, lastmeasstate) = oldtrajectory.findFirstLastMeasurement();
6871 std::unique_ptr<const TrackParameters> per(
nullptr);
6874 std::unique_ptr<const TrackParameters> prevpar(
6875 firstmeasstate->trackParameters() !=
nullptr ?
6876 firstmeasstate->trackParameters()->clone() :
6879 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.upstreamMaterialLayers();
6882 for (
int i = (
int)upstreamlayers.size() - 1;
i >= 0;
i--) {
6883 if (prevpar ==
nullptr) {
6890 if (
layer ==
nullptr) {
6891 layer = upstreamlayers[
i].second;
6894 DistanceSolution distsol =
layer->surfaceRepresentation().straightLineDistanceEstimate(
6895 prevpar->position(), prevpar->momentum().unit()
6899 if (distsol.numberOfSolutions() == 2) {
6904 if (distsol.first() * distsol.second() < 0 && !
first) {
6913 std::unique_ptr<const TrackParameters> layerpar(
6917 layer->surfaceRepresentation(),
6920 oldtrajectory.m_fieldprop,
6925 if (layerpar ==
nullptr) {
6929 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
6933 prevpar = std::move(layerpar);
6937 const Layer *startlayer = firstmeasstate->trackParameters()->associatedSurface().associatedLayer();
6939 if ((startlayer !=
nullptr) && (startlayer->layerMaterialProperties() !=
nullptr)) {
6940 double startfactor = startlayer->layerMaterialProperties()->alongPostFactor();
6941 const Surface & discsurf = startlayer->surfaceRepresentation();
6944 startfactor = startlayer->layerMaterialProperties()->oppositePostFactor();
6946 if (startfactor > 0.5) {
6947 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6948 firstmeasstate->trackParameters(), *startlayer,
oppositeMomentum, matEffects
6951 if (updatedpar !=
nullptr) {
6952 firstmeasstate->setTrackParameters(std::move(updatedpar));
6961 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
6963 if ((endlayer !=
nullptr) && (endlayer->layerMaterialProperties() !=
nullptr)) {
6964 double endfactor = endlayer->layerMaterialProperties()->alongPreFactor();
6965 const Surface & discsurf = endlayer->surfaceRepresentation();
6968 endfactor = endlayer->layerMaterialProperties()->oppositePreFactor();
6971 if (endfactor > 0.5) {
6972 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6973 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
6976 if (updatedpar !=
nullptr) {
6977 lastmeasstate->setTrackParameters(std::move(updatedpar));
6982 if (prevpar !=
nullptr) {
6989 oldtrajectory.m_fieldprop,
6994 if (per ==
nullptr) {
6995 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7000 }
else if (cache.m_acceleration && (firstmeasstate->trackParameters() !=
nullptr)) {
7002 *firstmeasstate->trackParameters(),
7008 per.reset(oldtrajectory.referenceParameters()->clone());
◆ myfit()
Definition at line 4626 of file GlobalChi2Fitter.cxx.
4634 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4638 if (!trajectory.m_straightline) {
4639 if (trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
4640 trajectory.m_straightline = !cache.m_field_cache.solenoidOn();
4641 }
else if ((trajectory.prefit() == 0) && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0) {
4642 trajectory.m_straightline = !cache.m_field_cache.toroidOn();
4644 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
4649 cache.m_lastiter = 0;
4653 if (trajectory.numberOfPerigeeParameters() == -1) {
4654 cache.incrementFitStatus(
S_FITS);
4655 if (trajectory.m_straightline) {
4656 trajectory.setNumberOfPerigeeParameters(4);
4658 trajectory.setNumberOfPerigeeParameters(5);
4662 if (trajectory.nDOF() < 0) {
4667 cache.m_phiweight.clear();
4668 cache.m_firstmeasurement.clear();
4669 cache.m_lastmeasurement.clear();
4672 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4676 if (matEffects ==
Trk::electron && trajectory.m_straightline) {
4682 trajectory.setMass(
mass);
4684 ATH_MSG_DEBUG(
"start param: " << param <<
" pos: " << param.position() <<
" pt: " << param.pT());
4686 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4688 if (!cache.m_acceleration && (per ==
nullptr)) {
4698 cache.m_acceleration &&
4699 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits() &&
4702 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4705 cache.m_fastmat =
false;
4710 !cache.m_acceleration ||
4711 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() != trajectory.numberOfHits()
4713 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4716 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4720 if (cache.m_acceleration && (trajectory.referenceParameters() ==
nullptr) && (per ==
nullptr)) {
4723 if (trajectory.numberOfScatterers() >= 2) {
4724 GXFTrackState *scatstate =
nullptr;
4725 GXFTrackState *scatstate2 =
nullptr;
4728 for (std::vector<std::unique_ptr<GXFTrackState>>::
iterator it =
4729 trajectory.trackStates().begin();
4730 it != trajectory.trackStates().end();
4734 scatindex == trajectory.numberOfScatterers() / 2 ||
4735 (**it).materialEffects()->deltaE() == 0
4737 scatstate2 = (*it).get();
4742 scatstate = (*it).get();
4749 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4750 throw std::logic_error(
"Invalid scatterer");
4753 vertex = .49 * (scatstate->position() + scatstate2->position());
4755 int nstates = (
int) trajectory.trackStates().size();
4757 trajectory.trackStates()[nstates / 2 - 1]->position() +
4758 trajectory.trackStates()[nstates / 2]->position()
4762 PerigeeSurface persurf(
vertex);
4763 std::unique_ptr<const TrackParameters> nearestpar;
4764 double mindist = 99999;
4765 std::vector < GXFTrackState * >mymatvec;
4767 for (
auto &
it : trajectory.trackStates()) {
4768 if ((*it).trackParameters() ==
nullptr) {
4773 .straightLineDistanceEstimate(
4774 (*it).trackParameters()->position(),
4775 (*it).trackParameters()->momentum().unit())
4779 (cache.m_caloEntrance ==
nullptr) ||
4780 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4784 (((*it).measurement() !=
nullptr) && insideid) || (
4785 ((*it).materialEffects() !=
nullptr) &&
4787 (*it).materialEffects()->deltaE() == 0 ||
4788 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4790 (*it).materialEffects()->deltaPhi() != 0
4794 double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4795 if (dist < mindist) {
4803 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4804 mymatvec.push_back(
it.get());
4808 if (nearestpar ==
nullptr) {
4812 for (
auto & state : mymatvec) {
4814 const Surface &matsurf = state->associatedSurface();
4815 DistanceSolution distsol = matsurf.straightLineDistanceEstimate(
4816 nearestpar->position(), nearestpar->momentum().unit());
4820 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4824 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4830 trajectory.m_fieldprop,
4834 if (tmppar ==
nullptr) {
4842 trajectory.m_fieldprop,
4846 if (tmppar ==
nullptr) {
4858 if (state->materialEffects()->sigmaDeltaE() > 0) {
4859 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4862 double de = std::abs(state->materialEffects()->deltaE());
4864 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4870 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4871 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4875 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4881 trajectory.m_fieldprop,
4886 if (tmpPars !=
nullptr) {
4887 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4892 double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4893 double toteloss = std::abs(trajectory.totalEnergyLoss());
4894 double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4898 per = per->associatedSurface().createUniqueTrackParameters(
4903 if (per ==
nullptr) {
4911 PerigeeSurface persurf2(per->position());
4912 per = persurf2.createUniqueTrackParameters(
4920 }
else if (per ==
nullptr) {
4924 if ((per ==
nullptr) && (trajectory.referenceParameters() ==
nullptr)) {
4932 if (trajectory.m_straightline && (per !=
nullptr)) {
4933 if (trajectory.numberOfPerigeeParameters() == -1) {
4934 trajectory.setNumberOfPerigeeParameters(4);
4938 per = per->associatedSurface().createUniqueTrackParameters(
4941 }
else if (trajectory.numberOfPerigeeParameters() == -1) {
4942 trajectory.setNumberOfPerigeeParameters(5);
4945 if (per !=
nullptr) {
4946 trajectory.setReferenceParameters(std::move(per));
4949 int nfitpar = trajectory.numberOfFitParameters();
4950 int nperpars = trajectory.numberOfPerigeeParameters();
4951 int nscat = trajectory.numberOfScatterers();
4952 int nbrem = trajectory.numberOfBrems();
4955 Eigen::MatrixXd a_inv;
4956 a.resize(nfitpar, nfitpar);
4961 derivPool.setZero();
4963 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
4964 if (state->materialEffects() !=
nullptr) {
4967 state->setDerivatives(derivPool);
4970 bool doderiv =
true;
4972 int tmpminiter = cache.m_miniter;
4975 cache.m_lastiter =
it;
4981 cache.m_miniter = tmpminiter;
4985 if (!trajectory.converged()) {
4986 cache.m_fittercode =
4996 cache.m_miniter = tmpminiter;
5000 int nhits = trajectory.numberOfHits();
5001 int ntrthits = trajectory.numberOfTRTHits();
5002 int nsihits = trajectory.numberOfSiliconHits();
5003 double redchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
5004 double prevredchi2 = (trajectory.nDOF() > 0) ? trajectory.prevchi2() / trajectory.nDOF() : 0;
5013 (redchi2 < prevredchi2 &&
5014 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
5015 nsihits + ntrthits == nhits
5020 if (
it != 1 || nsihits != 0 || trajectory.nDOF() <= 0 || trajectory.chi2() / trajectory.nDOF() <= 3) {
5025 cache.m_miniter = tmpminiter;
5032 int ntrtprechits = trajectory.numberOfTRTPrecHits();
5033 int ntrttubehits = trajectory.numberOfTRTTubeHits();
5035 if (ntrtprechits+ntrttubehits) {
5036 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
5038 if (phf<m_minphfcut && it>=3) {
5039 if ((ntrtprechits+ntrttubehits)>=15) {
5044 <<
" | nTRTPrecHits = " << ntrtprechits
5045 <<
" | nTRTTubeHits = " << ntrttubehits
5046 <<
" | nOutliers = "
5047 << trajectory.numberOfOutliers());
5049 if (!trajectory.converged()) {
5055 cache.m_miniter = tmpminiter;
5064 cache.m_miniter = tmpminiter;
5066 if (trajectory.prefit() == 0) {
5070 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
5071 if (lltOfW.info() == Eigen::Success) {
5075 int ncols =
a.cols();
5076 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
5077 a_inv = lltOfW.solve(weightInvAMG);
5086 GXFTrajectory *finaltrajectory = &trajectory;
5088 (runOutlier || cache.m_sirecal) &&
5089 trajectory.numberOfSiliconHits() == trajectory.numberOfHits()
5097 if (finaltrajectory != &trajectory) {
5098 delete finaltrajectory;
5107 if (!cache.m_acceleration && (finaltrajectory->prefit() == 0)) {
5108 if (nperpars == 5) {
5109 for (
int i = 0;
i <
a.cols();
i++) {
5110 a_inv(4,
i) *= .001;
5111 a_inv(
i, 4) *= .001;
5115 int scatterPos = nperpars + 2 * nscat;
5116 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
5117 for (
int i = 0;
i <
a.cols();
i++) {
5118 a_inv(scatterPos,
i) *= .001;
5119 a_inv(
i, scatterPos) *= .001;
5125 int nperparams = finaltrajectory->numberOfPerigeeParameters();
5126 for (
int i = 0;
i < nperparams;
i++) {
5127 for (
int j = 0; j < nperparams; j++) {
5128 (errmat) (j,
i) = a_inv(j,
i);
5132 if (trajectory.m_straightline) {
5133 (errmat) (4, 4) = 1
e-20;
5136 const AmgVector(5) & perpars = finaltrajectory->referenceParameters()->parameters();
5137 std::unique_ptr<const TrackParameters> measper(
5138 finaltrajectory->referenceParameters()->associatedSurface().createUniqueTrackParameters(
5139 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5143 finaltrajectory->setReferenceParameters(std::move(measper));
5145 cache.m_fullcovmat = a_inv;
5149 std::unique_ptr<Track>
track =
nullptr;
5151 if (finaltrajectory->prefit() > 0) {
5152 if (finaltrajectory != &trajectory) {
5153 delete finaltrajectory;
5158 if (finaltrajectory->numberOfOutliers() <=
m_maxoutliers || !runOutlier) {
5165 double cut = (finaltrajectory->numberOfSiliconHits() ==
5166 finaltrajectory->numberOfHits())
5172 (
track !=
nullptr) && (
5173 track->fitQuality()->numberDoF() != 0 &&
5174 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5177 track.reset(
nullptr);
5181 if (
track ==
nullptr) {
5185 if (finaltrajectory != &trajectory) {
5186 delete finaltrajectory;
5189 return track.release();
◆ myfit_helper()
◆ numericalDerivatives()
Definition at line 8123 of file GlobalChi2Fitter.cxx.
8130 ParamDefsAccessor paraccessor;
8138 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8141 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8148 const Surface & previousSurface = tmpprevpar->associatedSurface();
8152 for (
int i = 0;
i < 5;
i++) {
8155 if (thisdiscsurf &&
i == 1) {
8159 vecpluseps[paraccessor.pardef[
i]] += eps[
i];
8160 vecminuseps[paraccessor.pardef[
i]] -= eps[
i];
8161 if (thiscylsurf &&
i == 0) {
8162 if (vecpluseps[0] / previousSurface.bounds().r() >
M_PI) {
8163 vecpluseps[0] -= 2 *
M_PI * previousSurface.bounds().r();
8165 if (vecminuseps[0] / previousSurface.bounds().r() < -
M_PI) {
8166 vecminuseps[0] += 2 *
M_PI * previousSurface.bounds().r();
8169 if (thisdiscsurf &&
i == 1) {
8170 if (vecpluseps[
i] >
M_PI) {
8171 vecpluseps[
i] -= 2 *
M_PI;
8173 if (vecminuseps[
i] < -
M_PI) {
8174 vecminuseps[
i] += 2 *
M_PI;
8180 std::unique_ptr<const TrackParameters> parpluseps(
8181 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8190 std::unique_ptr<const TrackParameters> parminuseps(
8191 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8201 std::unique_ptr<const TrackParameters> newparpluseps(
8212 std::unique_ptr<const TrackParameters> newparminuseps(
8227 if (newparpluseps ==
nullptr) {
8239 if (newparminuseps ==
nullptr) {
8251 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8255 for (
int j = 0; j < 5; j++) {
8256 double diff = newparpluseps->parameters()[paraccessor.pardef[j]] -
8257 newparminuseps->parameters()[paraccessor.pardef[j]];
8258 if (cylsurf && j == 0) {
8259 double length = 2 *
M_PI * surf.bounds().r();
8268 if (discsurf && j == 1) {
8269 if (std::abs(std::abs(
diff) - 2 *
M_PI) < std::abs(
diff)) {
8278 (*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 914 of file GlobalChi2Fitter.h.
919 if (!handle.isValid()) {
922 return handle.cptr();
◆ runIteration()
Definition at line 5660 of file GlobalChi2Fitter.cxx.
5671 int nfitpars = trajectory.numberOfFitParameters();
5672 int nperpars = trajectory.numberOfPerigeeParameters();
5673 int scatpars = 2 * trajectory.numberOfScatterers();
5674 int nupstreamstates = trajectory.numberOfUpstreamStates();
5675 int nbrem = trajectory.numberOfBrems();
5676 double oldchi2 = trajectory.chi2();
5677 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5678 int nsihits = trajectory.numberOfSiliconHits();
5679 int ntrthits = trajectory.numberOfTRTHits();
5680 int nhits = trajectory.numberOfHits();
5682 if (cache.m_phiweight.empty()) {
5683 cache.m_phiweight.assign(trajectory.trackStates().size(), 1);
5696 double newredchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
5698 ATH_MSG_DEBUG(
"old chi2: " << oldchi2 <<
"/" << trajectory.nDOF() <<
5699 "=" << oldredchi2 <<
" new chi2: " << trajectory.chi2() <<
"/" <<
5700 trajectory.nDOF() <<
"=" << newredchi2);
5702 if (trajectory.prefit() > 0 && trajectory.converged()) {
5708 std::vector < std::pair < double, double >>&scatsigmas = trajectory.scatteringSigmas();
5710 int nmeas = (
int)
res.size();
5712 const Amg::MatrixX & weight_deriv = trajectory.weightedResidualDerivatives();
5719 if (cache.m_firstmeasurement.empty()) {
5720 cache.m_firstmeasurement.resize(nfitpars);
5721 cache.m_lastmeasurement.resize(nfitpars);
5722 for (
int i = 0;
i < nperpars;
i++) {
5723 cache.m_firstmeasurement[
i] = 0;
5724 cache.m_lastmeasurement[
i] = nmeas - nbrem;
5729 for (
int i = 0;
i < (
int) trajectory.trackStates().size();
i++) {
5730 std::unique_ptr<GXFTrackState> & state = trajectory.trackStates()[
i];
5731 GXFMaterialEffects *meff = state->materialEffects();
5732 if (meff ==
nullptr) {
5733 measno += state->numberOfMeasuredParameters();
5735 if (meff !=
nullptr) {
5736 if (meff->sigmaDeltaTheta() != 0
5737 && ((trajectory.prefit() == 0) || meff->deltaE() == 0)) {
5738 int scatterPos = nperpars + 2 * scatno;
5739 if (
i < nupstreamstates) {
5740 cache.m_lastmeasurement[scatterPos] =
5741 cache.m_lastmeasurement[scatterPos + 1] = measno;
5742 cache.m_firstmeasurement[scatterPos] =
5743 cache.m_firstmeasurement[scatterPos + 1] = 0;
5745 cache.m_lastmeasurement[scatterPos] =
5746 cache.m_lastmeasurement[scatterPos + 1] = nmeas - nbrem;
5747 cache.m_firstmeasurement[scatterPos] =
5748 cache.m_firstmeasurement[scatterPos + 1] = measno;
5752 if (meff->sigmaDeltaE() > 0) {
5753 if (
i < nupstreamstates) {
5754 cache.m_firstmeasurement[nperpars + scatpars + bremno] = 0;
5755 cache.m_lastmeasurement[nperpars + scatpars + bremno] = measno;
5757 cache.m_firstmeasurement[nperpars + scatpars + bremno] = measno;
5758 cache.m_lastmeasurement[nperpars + scatpars + bremno] =
5768 if (
a.cols() != nfitpars) {
5772 for (
int k = 0;
k < nfitpars;
k++) {
5774 int maxmeas = nmeas - nbrem;
5775 maxmeas = cache.m_lastmeasurement[
k];
5776 minmeas = cache.m_firstmeasurement[
k];
5778 for (measno = minmeas; measno < maxmeas; measno++) {
5780 res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5784 if (
k == 4 ||
k >= nperpars + scatpars) {
5785 for (measno = nmeas - nbrem; measno < nmeas; measno++) {
5786 b[
k] +=
res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5791 for (
int l =
k;
l < nfitpars;
l++) {
5793 std::min(cache.m_lastmeasurement[
k], cache.m_lastmeasurement[
l]);
5796 cache.m_firstmeasurement[
l]);
5798 for (measno = minmeas; measno < maxmeas; measno++) {
5799 tmp += weight_deriv(measno,
k) * weight_deriv(measno,
l);
5801 a.fillSymmetric(
l,
k,
tmp);
5809 for (
int k = nperpars;
k < nperpars + scatpars;
k += 2) {
5810 a(
k,
k) += 1. / (scatsigmas[scatno].first * scatsigmas[scatno].first);
5811 a(
k + 1,
k + 1) += 1. / (scatsigmas[scatno].second * scatsigmas[scatno].second);
5815 for (
int measno = nmeas - nbrem; measno < nmeas; measno++) {
5816 for (
int k = 4;
k < nfitpars;
k++) {
5818 k = nperpars + scatpars;
5821 for (
int l =
k;
l < nfitpars;
l++) {
5823 l = nperpars + scatpars;
5825 double tmp =
a(
l,
k) + weight_deriv(measno,
k) * weight_deriv(measno,
l);
5826 a.fillSymmetric(
l,
k,
tmp);
5832 unsigned int scatno = 0;
5833 bool weightchanged =
false;
5835 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5836 GXFMaterialEffects *meff = thisstate->materialEffects();
5838 if (meff !=
nullptr) {
5839 const PlaneSurface *plsurf =
nullptr;
5842 plsurf =
static_cast < const PlaneSurface *
>(&thisstate->associatedSurface());
5843 if (meff->deltaE() == 0 || ((trajectory.prefit() == 0) && (plsurf !=
nullptr))) {
5844 weightchanged =
true;
5846 if (
a.cols() != nfitpars) {
5850 int scatNoIndex = 2 * scatno + nperpars;
5852 if (trajectory.prefit() == 0) {
5853 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5854 if (scatno >= cache.m_phiweight.size()) {
5856 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5857 throw std::range_error(
message.str());
5861 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5865 cache.m_phiweight[scatno] = 1.00000001;
5866 }
else if (
it == 1) {
5867 cache.m_phiweight[scatno] = 1.0000001;
5868 }
else if (
it <= 3) {
5869 cache.m_phiweight[scatno] = 1.0001;
5870 }
else if (
it <= 6) {
5871 cache.m_phiweight[scatno] = 1.01;
5873 cache.m_phiweight[scatno] = 1.1;
5876 a(scatNoIndex, scatNoIndex) *= cache.m_phiweight[scatno];
5880 else if (trajectory.prefit() >= 2) {
5881 if (newredchi2 > oldredchi2 - 1 && newredchi2 < oldredchi2) {
5882 a(scatNoIndex, scatNoIndex) *= 1.0001;
5883 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5884 }
else if (newredchi2 > oldredchi2 - 25 && newredchi2 < oldredchi2) {
5885 a(scatNoIndex, scatNoIndex) *= 1.001;
5886 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5888 a(scatNoIndex, scatNoIndex) *= 1.1;
5889 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.001;
5895 thisstate->materialEffects()->sigmaDeltaPhi() != 0 &&
5896 ((trajectory.prefit() == 0) || thisstate->materialEffects()->deltaE() == 0)
5904 (trajectory.prefit() == 2) &&
5906 trajectory.numberOfBrems() > 0 &&
5907 (newredchi2 < oldredchi2 - 25 || newredchi2 > oldredchi2)
5912 if (doderiv || weightchanged) {
5916 if (trajectory.converged()) {
5917 if ((trajectory.prefit() == 0) && nsihits + ntrthits != nhits) {
5918 unsigned int scatno = 0;
5920 if (
a.cols() != nfitpars) {
5924 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5925 if ((thisstate->materialEffects() !=
nullptr) && thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5926 if (scatno >= cache.m_phiweight.size()) {
5928 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5929 throw std::range_error(
message.str());
5932 const PlaneSurface *plsurf =
nullptr;
5935 plsurf =
static_cast<const PlaneSurface *
>(&thisstate->associatedSurface());
5937 if (thisstate->materialEffects()->deltaE() == 0 || (plsurf !=
nullptr)) {
5938 int scatNoIndex = 2 * scatno + nperpars;
5939 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5940 cache.m_phiweight[scatno] = 1;
5943 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5956 (newredchi2 < 2 || (newredchi2 < oldredchi2 && newredchi2 > oldredchi2 - .5)) &&
5957 (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 6320 of file GlobalChi2Fitter.cxx.
6329 bool trackok =
false;
6330 GXFTrajectory *oldtrajectory = &trajectory;
6331 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6332 GXFTrajectory *newtrajectory =
nullptr;
6333 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6339 while (!trackok && oldtrajectory->nDOF() > 0) {
6341 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->trackStates();
6344 Amg::MatrixX & weightderiv = oldtrajectory->weightedResidualDerivatives();
6345 int nfitpars = oldtrajectory->numberOfFitParameters();
6346 int nhits = oldtrajectory->numberOfHits();
6347 int nsihits = oldtrajectory->numberOfSiliconHits();
6349 if (nhits != nsihits) {
6353 double maxsipull = -1;
6355 int hitno_maxsipull = -1;
6356 int measno_maxsipull = -1;
6357 int stateno_maxsipull = 0;
6358 GXFTrackState *state_maxsipull =
nullptr;
6363 int noutl = oldtrajectory->numberOfOutliers();
6369 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6370 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6376 double *
errors = state->measurementErrors();
6378 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6379 double sinstereo = state->sinStereo();
6380 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6381 double weight1 = -1;
6383 if (hitcov(0, 0) > trackcov(0, 0)) {
6384 if (sinstereo == 0) {
6388 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6389 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6400 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6403 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6406 sipull1 =
std::max(sipull1, sipull2);
6408 if (sipull1 > maxsipull) {
6409 maxsipull = sipull1;
6410 measno_maxsipull = measno;
6411 state_maxsipull = state.get();
6412 stateno_maxsipull = stateno;
6413 hitno_maxsipull = hitno;
6424 measno += state->numberOfMeasuredParameters();
6428 double maxpull = maxsipull;
6430 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6431 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6440 oldtrajectory->chi2() / oldtrajectory->nDOF() > .25 *
m_chi2cut
6442 state_maxsipull = oldtrajectory->trackStates()[stateno_maxsipull].get();
6443 const PrepRawData *prd{};
6445 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6446 prd = rot->prepRawData();
6448 std::unique_ptr < const RIO_OnTrack > broadrot;
6449 double *olderror = state_maxsipull->measurementErrors();
6451 const TrackParameters *trackpar_maxsipull = state_maxsipull->trackParameters();
6453 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6454 std::unique_ptr<const TrackParameters> trackparForCorrect(
6455 trackpar_maxsipull->associatedSurface().createUniqueTrackParameters(
6461 state_maxsipull->hasTrackCovariance()
6463 state_maxsipull->trackCovariance())
6467 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6468 double newpull = -1;
6469 double newpull1 = -1;
6470 double newpull2 = -1;
6471 double newres1 = -1;
6472 double newres2 = -1;
6473 double newsinstereo = 0;
6477 !state_maxsipull->isRecalibrated() &&
6479 oldtrajectory->chi2() / trajectory.nDOF() > .3 *
m_chi2cut &&
6487 newerror[0] = std::sqrt(
covmat(0, 0));
6489 if (state_maxsipull->sinStereo() != 0) {
6507 newerror[0] = std::sqrt(
v0);
6510 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6512 if (cosstereo != 1.) {
6514 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6515 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6518 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6521 if (newerror[0] == 0.0) {
6522 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6526 newpull1 = std::abs(newres1 / newerror[0]);
6530 newerror[1] = std::sqrt(
covmat(1, 1));
6531 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6532 newpull2 = std::abs(newres2 / newerror[1]);
6535 newpull =
std::max(newpull1, newpull2);
6541 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6543 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6544 throw std::runtime_error(
6545 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6550 newtrajectory = oldtrajectory;
6552 if (
a.cols() != nfitpars) {
6556 double oldres1 =
res[measno_maxsipull];
6557 res[measno_maxsipull] = newres1;
6558 err[measno_maxsipull] = newerror[0];
6560 for (
int i = 0;
i < nfitpars;
i++) {
6561 if (weightderiv(measno_maxsipull,
i) == 0) {
6565 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6567 for (
int j =
i; j < nfitpars; j++) {
6571 weightderiv(measno_maxsipull,
i) *
6572 weightderiv(measno_maxsipull, j) *
6573 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6577 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6581 double oldres2 =
res[measno_maxsipull + 1];
6582 res[measno_maxsipull + 1] = newres2;
6583 err[measno_maxsipull + 1] = newerror[1];
6585 for (
int i = 0;
i < nfitpars;
i++) {
6586 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6590 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6592 for (
int j =
i; j < nfitpars; j++) {
6596 weightderiv(measno_maxsipull + 1,
i) *
6597 weightderiv(measno_maxsipull + 1, j) *
6598 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6603 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6608 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6609 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6610 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6611 olderror[1] <<
" newerror_1=" << newerror[1]
6614 state_maxsipull->setMeasurement(std::move(broadrot));
6615 state_maxsipull->setSinStereo(newsinstereo);
6616 state_maxsipull->setMeasurementErrors(newerror);
6620 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6621 (oldtrajectory->chi2() / oldtrajectory->nDOF() > .3 *
m_chi2cut || noutl > 1)
6625 (oldtrajectory->nDOF() > 1 || hittype_maxsipull ==
TrackState::SCT) &&
6630 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6631 measno_maxsipull <<
" pull=" << maxsipull
6638 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6639 newtrajectory = cleanup_newtrajectory.get();
6641 if (newa.cols() != nfitpars) {
6646 Amg::MatrixX & newweightderiv = newtrajectory->weightedResidualDerivatives();
6647 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6648 throw std::runtime_error(
6649 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6653 double oldres1 =
res[measno_maxsipull];
6654 newres[measno_maxsipull] = 0;
6656 for (
int i = 0;
i < nfitpars;
i++) {
6657 if (weightderiv(measno_maxsipull,
i) == 0) {
6661 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6663 for (
int j =
i; j < nfitpars; j++) {
6667 weightderiv(measno_maxsipull,
i) *
6668 weightderiv(measno_maxsipull, j)
6672 newweightderiv(measno_maxsipull,
i) = 0;
6676 double oldres2 =
res[measno_maxsipull + 1];
6677 newres[measno_maxsipull + 1] = 0;
6679 for (
int i = 0;
i < nfitpars;
i++) {
6680 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6684 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6686 for (
int j =
i; j < nfitpars; j++) {
6687 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6694 weightderiv(measno_maxsipull + 1,
i) *
6695 weightderiv(measno_maxsipull + 1, j)
6699 newweightderiv(measno_maxsipull + 1,
i) = 0;
6703 newtrajectory->setOutlier(stateno_maxsipull);
6709 newtrajectory->setConverged(
false);
6725 if (!newtrajectory->converged()) {
6727 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6734 if (!newtrajectory->converged()) {
6743 double oldchi2 = oldtrajectory->chi2() / oldtrajectory->nDOF();
6744 double newchi2 = (newtrajectory->nDOF() > 0) ? newtrajectory->chi2() / newtrajectory->nDOF() : 0;
6747 if (newtrajectory->nDOF() != oldtrajectory->nDOF() && maxsipull > cut2) {
6748 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6750 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6755 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6756 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6764 (void)cleanup_oldtrajectory.release();
6765 return oldtrajectory;
6767 if (oldtrajectory != newtrajectory) {
6768 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6769 oldtrajectory = newtrajectory;
6776 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6777 if (lltOfW.info() == Eigen::Success) {
6781 int ncols =
a.cols();
6782 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6783 fullcov = lltOfW.solve(weightInvAMG);
6801 oldtrajectory->nDOF() > 0 &&
6802 oldtrajectory->chi2() / oldtrajectory->nDOF() >
m_chi2cut &&
6810 (void)cleanup_oldtrajectory.release();
6811 return oldtrajectory;
◆ runTrackCleanerTRT()
Definition at line 6150 of file GlobalChi2Fitter.cxx.
6163 if (
it == 1 && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
6167 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6170 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6171 int nfitpars = trajectory.numberOfFitParameters();
6173 if (
a.cols() != nfitpars) {
6177 int nperpars = trajectory.numberOfPerigeeParameters();
6178 int nscats = trajectory.numberOfScatterers();
6181 bool outlierremoved =
false;
6182 bool hitrecalibrated =
false;
6184 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6185 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6193 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6197 trajectory.setOutlier(stateno);
6198 outlierremoved =
true;
6200 double *
errors = state->measurementErrors();
6201 double olderror =
errors[0];
6203 trajectory.updateTRTHitCount(stateno, olderror);
6205 for (
int i = 0;
i < nfitpars;
i++) {
6206 if (weightderiv(measno,
i) == 0) {
6210 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6212 for (
int j =
i; j < nfitpars; j++) {
6215 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6218 weightderiv(measno,
i) = 0;
6222 }
else if (trtrecal) {
6223 double *
errors = state->measurementErrors();
6224 double olderror =
errors[0];
6226 const auto *
const thisMeasurement{state->measurement()};
6232 if (oldrot->prepRawData() !=
nullptr) {
6233 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6235 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6237 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6238 double distance = std::abs(std::abs(trackradius) - dcradius);
6240 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6241 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6242 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6243 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6246 if (newrot !=
nullptr) {
6248 hitrecalibrated =
true;
6252 if ((measno < 0) or (measno >= (
int)
res.size())) {
6253 throw std::runtime_error(
6254 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6258 double oldres =
res[measno];
6259 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6261 state->setMeasurement(std::move(newrot));
6263 trajectory.updateTRTHitCount(stateno, olderror);
6265 for (
int i = 0;
i < nfitpars;
i++) {
6266 if (weightderiv(measno,
i) == 0) {
6270 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6272 for (
int j =
i; j < nfitpars; j++) {
6276 !cache.m_phiweight.empty() &&
6279 i < nperpars + 2 * nscats &&
6280 (
i - nperpars) % 2 == 0
6282 weight = cache.m_phiweight[(
i - nperpars) / 2];
6287 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6290 weightderiv(measno,
i) *= olderror / newerror;
6293 res[measno] = newres;
6294 err[measno] = newerror;
6303 measno += state->numberOfMeasuredParameters();
6307 if (trajectory.nDOF() < 0) {
6312 if (outlierremoved || hitrecalibrated) {
6314 trajectory.setConverged(
false);
6316 cache.m_miniter =
it + 2;
◆ setMinIterations()
void Trk::GlobalChi2Fitter::setMinIterations |
( |
int |
| ) |
|
|
privatevirtual |
Definition at line 8289 of file GlobalChi2Fitter.cxx.
8291 (
"Configure the minimum number of Iterations via jobOptions");
◆ throwFailedToGetTrackingGeomtry()
void Trk::GlobalChi2Fitter::throwFailedToGetTrackingGeomtry |
( |
| ) |
const |
|
private |
◆ trackingGeometry()
Definition at line 907 of file GlobalChi2Fitter.h.
910 if (!cache.m_trackingGeometry)
912 return cache.m_trackingGeometry;
◆ updateEnergyLoss()
Definition at line 7787 of file GlobalChi2Fitter.cxx.
7800 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7804 double newqoverp = 0;
7806 if (meff.sigmaDeltaE() <= 0) {
7811 double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.deltaE()) * std::sqrt(
mass *
mass + oldp * oldp) + meff.deltaE() * meff.deltaE();
7818 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7824 return surf.createUniqueTrackParameters(
7825 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
◆ updateFitParameters()
Definition at line 5965 of file GlobalChi2Fitter.cxx.
5973 double d0 = refpar->parameters()[
Trk::d0];
5974 double z0 = refpar->parameters()[
Trk::z0];
5977 double qoverp = refpar->parameters()[
Trk::qOverP];
5978 int nscat = trajectory.numberOfScatterers();
5979 int nbrem = trajectory.numberOfBrems();
5980 int nperparams = trajectory.numberOfPerigeeParameters();
5982 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
5985 if (llt.info() == Eigen::Success) {
5991 if (trajectory.numberOfPerigeeParameters() > 0) {
5996 qoverp = (trajectory.m_straightline) ? 0 : .001 *
result[4] + qoverp;
6005 std::vector < std::pair < double, double >>&scatangles = trajectory.scatteringAngles();
6006 std::vector < double >&delta_ps = trajectory.brems();
6008 for (
int i = 0;
i < nscat;
i++) {
6009 scatangles[
i].first +=
result[2 *
i + nperparams];
6010 scatangles[
i].second +=
result[2 *
i + nperparams + 1];
6013 for (
int i = 0;
i < nbrem;
i++) {
6014 delta_ps[
i] +=
result[nperparams + 2 * nscat +
i];
6017 std::unique_ptr<const TrackParameters> newper(
6018 trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
6023 trajectory.setReferenceParameters(std::move(newper));
6024 trajectory.setScatteringAngles(scatangles);
6025 trajectory.setBrems(delta_ps);
◆ updatePixelROTs()
Update the Pixel ROT using the current trajectory/local track parameters.
Definition at line 6030 of file GlobalChi2Fitter.cxx.
6036 if ( trajectory.numberOfSiliconHits() == 0) {
6045 if (!splitProbContainer.isValid()) {
6049 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6052 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6053 int nfitpars = trajectory.numberOfFitParameters();
6056 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6061 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6064 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6074 const PrepRawData *prd{};
6076 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6077 prd = rot->prepRawData();
6087 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6088 if (!splitProb.isSplit()) {
6089 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6093 std::unique_ptr < const RIO_OnTrack > newrot;
6094 double *olderror = state->measurementErrors();
6097 double newerror[5] = {-1,-1,-1,-1,-1};
6098 double newres[2] = {-1,-1};
6100 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6107 newerror[0] = std::sqrt(
covmat(0, 0));
6108 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6109 newerror[1] = std::sqrt(
covmat(1, 1));
6110 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6112 if (
a.cols() != nfitpars) {
6117 for(
int k =0;
k<2;
k++ ){
6118 double oldres =
res[measno+
k];
6119 res[measno+
k] = newres[
k];
6120 err[measno+
k] = newerror[
k];
6122 for (
int i = 0;
i < nfitpars;
i++) {
6123 if (weightderiv(measno+
k,
i) == 0) {
6127 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6129 for (
int j =
i; j < nfitpars; j++) {
6133 weightderiv(measno+
k,
i) *
6134 weightderiv(measno+
k, j) *
6135 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6139 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6143 state->setMeasurement(std::move(newrot));
6144 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 932 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 925 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
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
void runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
void addMaterialUpdateTrajectory(Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * >> &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const
Given layer information, probe those layers for scatterers and add them to a track.
std::vector< size_t > vec
bool is_trt(Identifier id) const
std::vector< std::unique_ptr< TrackParameters > > holesearchExtrapolation(const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const
Helper method which performs an extrapolation with additional logic for hole search.
#define ATH_MSG_VERBOSE(x)
ToolHandle< IRIO_OnTrackCreator > m_broadROTcreator
base class to integrate material effects on Trk::Track in a flexible way.
std::unique_ptr< const TrackParameters > makeTrackFindPerigeeParameters(const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const
bool empty() const
Test if the key is blank.
Gaudi::Property< int > m_maxoutliers
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ numberOfSCTHoles
number of Holes in both sides of a SCT module
bool is_valid() const
Check if id is in a valid state.
SG::ReadCondHandleKey< TrackingGeometry > m_trackingGeometryReadKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_field_cache_key
ToolHandle< IBoundaryCheckTool > m_boundaryCheckTool
@ u
Enums for curvilinear frames.
FitterStatusCode calculateTrackParameters(const EventContext &ctx, GXFTrajectory &, bool) const
void stable_sort(std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, std::reverse_iterator< DataModel_detail::iterator< DVL > > end, Compare comp)
Specialization of stable_sort for DataVector/List.
Gaudi::Property< int > m_fixbrem
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 ...
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
const LayerArray * confinedLayers() const
Return the subLayer array.
Gaudi::Property< int > m_maxit
@ StraightTrack
A straight track.
AmgSymMatrix(5) &GXFTrackState
@ SlimmedTrack
A slimmed track.
represents the full description of deflection and e-loss of a track in material.
ToolHandle< IResidualPullCalculator > m_residualPullCalculator
std::vector< std::reference_wrapper< GXFTrackState > > holeSearchStates(GXFTrajectory &trajectory) const
Extracts a collection of track states which are important for hole search.
static void addMaterialGetLayers(Cache &cache, std::vector< std::pair< const Layer *, const Layer * >> &layers, std::vector< std::pair< const Layer *, const Layer * >> &uplayers, const std::vector< std::unique_ptr< GXFTrackState >> &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat)
Collect all possible layers that a given track could have passed through.
static void makeTrackFillDerivativeMatrix(Cache &, GXFTrajectory &)
ToolHandle< Trk::ITrkMaterialProviderTool > m_caloMaterialProvider
double getDistance(const xAOD::Vertex *vtx1, const xAOD::Vertex *vtx2)
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 updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const
Update the Pixel ROT using the current trajectory/local track parameters.
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.
ParametersBase< TrackParametersDim, Charged > TrackParameters
bool is_indet(Identifier id) const
double charge(const T &p)
Gaudi::Property< double > m_chi2cut
StatusCode initialize(bool used=true)
Eigen::Matrix< double, 3, 1 > Vector3D
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
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 point B' along the line B that's closest to a second line A.
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
#define ATH_MSG_WARNING(x)
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
def delta_phi(phi1, phi2)
Gaudi::Property< int > m_maxitPixelROT
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
Track * fitIm(const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const
std::unique_ptr< T > unique_clone(const T *v)
Gaudi::Property< double > m_scalefactor
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
ToolHandle< IUpdator > m_updator
void addIDMaterialFast(const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final
@ PseudoMeasurementOnTrack
bool is_stgc(Identifier id) const
Gaudi::Property< bool > m_calomat
Gaudi::Property< bool > m_rejectLargeNScat
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
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
@ Unknown
Track fitter not defined.