|
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 7822 of file GlobalChi2Fitter.cxx.
7823 int nstatesupstream = trajectory.numberOfUpstreamStates();
7824 int nscatupstream = trajectory.numberOfUpstreamScatterers();
7825 int nbremupstream = trajectory.numberOfUpstreamBrems();
7826 int nscats = trajectory.numberOfScatterers();
7827 int nperpars = trajectory.numberOfPerigeeParameters();
7828 int nfitpars = trajectory.numberOfFitParameters();
7830 using Matrix55 = Eigen::Matrix<double, 5, 5>;
7832 Matrix55 initialjac;
7833 initialjac.setZero();
7834 initialjac(4, 4) = 1;
7836 Matrix55 jacvertex(initialjac);
7838 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.numberOfScatterers(), initialjac);
7839 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.numberOfBrems(), initialjac);
7841 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7842 GXFTrackState *prevstate =
nullptr, *state =
nullptr;
7844 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
7846 for (
bool forward : {
false,
true}) {
7848 hit_begin = nstatesupstream;
7850 scatno = nscatupstream;
7851 bremno = nbremupstream;
7853 hit_begin = nstatesupstream - 1;
7855 scatno = trajectory.numberOfUpstreamScatterers() - 1;
7856 bremno = trajectory.numberOfUpstreamBrems() - 1;
7860 int hitno = hit_begin;
7861 forward ? (hitno < hit_end) : (hitno >= hit_end);
7862 hitno += (forward ? 1 : -1)
7865 state =
states[hitno].get();
7869 if (fillderivmat && state->derivatives().cols() != nfitpars) {
7870 state->derivatives().resize(5, nfitpars);
7871 state->derivatives().setZero();
7874 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
7876 if (hitno == (forward ? hit_end - 1 : 0)) {
7877 if (!fillderivmat) {
7885 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
7887 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
7888 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
7889 jacvertex(4, 4) = jac(4, 4);
7891 int jmin = 0, jmax = 0, jcnt = 0;
7892 int lp_bgn = 0, lp_end = 0;
7896 jcnt = jmax - jmin + 1;
7898 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
7901 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7903 i == scatno + (forward ? -1 : 1) &&
7904 prevstate !=
nullptr &&
7906 (!trajectory.prefit() || prevstate->materialEffects()->deltaE() == 0)
7908 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7909 jacscat[
i](4, 4) = jac(4, 4);
7911 calculateJac(jac, jacscat[
i], jmin, jmax);
7915 Eigen::MatrixXd & derivmat = state->derivatives();
7916 int scatterPos = nperpars + 2 *
i;
7918 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
7924 jcnt = jmax - jmin + 1;
7926 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
7929 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
7931 i == bremno + (forward ? -1 : 1) &&
7933 prevstate->materialEffects() &&
7934 prevstate->materialEffects()->sigmaDeltaE() > 0
7936 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
7937 jacbrem[
i](4, 4) = jac(4, 4);
7939 calculateJac(jac, jacbrem[
i], jmin, jmax);
7943 Eigen::MatrixXd & derivmat = state->derivatives();
7944 int scatterPos = nperpars + 2 * nscats +
i;
7946 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
7950 calculateJac(jac, jacvertex, 0, 4);
7954 Eigen::MatrixXd & derivmat = state->derivatives();
7955 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
7957 if (nperpars == 5) {
7958 derivmat.col(4).segment(0, 4) *= .001;
7959 derivmat(4, 4) = .001 * jacvertex(4, 4);
7965 (!trajectory.prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
7967 scatno += (forward ? 1 : -1);
7971 states[hitno]->materialEffects() &&
7972 states[hitno]->materialEffects()->sigmaDeltaE() > 0
7974 bremno += (forward ? 1 : -1);
7977 prevstate =
states[hitno].get();
◆ calculateTrackErrors()
Definition at line 7984 of file GlobalChi2Fitter.cxx.
7992 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7993 int nstatesupstream = trajectory.numberOfUpstreamStates();
7995 GXFTrackState *prevstate =
nullptr;
7996 int i = nstatesupstream;
7997 for (
int j = 0; j < (
int)
states.size(); j++) {
7998 if (j < nstatesupstream) {
8005 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8006 if (stateno == 0 || stateno == nstatesupstream) {
8007 prevstate =
nullptr;
8010 std::unique_ptr<GXFTrackState> & state =
states[
index];
8011 if (state->materialEffects() !=
nullptr) {
8012 prevstate = state.get();
8016 if (!state->hasTrackCovariance()) {
8017 state->zeroTrackCovariance();
8019 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8021 if ((prevstate !=
nullptr) &&
8025 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8028 trackerrmat = jac * prevcov * jac.transpose();
8032 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8036 const MeasurementBase *measurement = state->measurement();
8037 const Amg::MatrixX & meascov = measurement->localCovariance();
8039 ParamDefsAccessor paraccessor;
8043 bool errorok =
true;
8044 for (
int i = 0;
i < 5;
i++) {
8045 if (measurement->localParameters().contains(paraccessor.pardef[
i])) {
8047 && trackerrmat(
i,
i) > meascov(j, j)) {
8049 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8050 trackerrmat(
i,
i) = meascov(j, j);
8051 for (
int k = 0;
k < 5;
k++) {
8061 for (
int i = 0;
i < 5;
i++) {
8065 for (
int j = 0; j < 5; j++) {
8072 if (trajectory.m_straightline) {
8073 trackerrmat(4, 4) = 1
e-20;
8077 state->trackParameters();
8079 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8081 if (state->hasTrackCovariance()) {
8082 trkerrmat = (state->trackCovariance());
8084 trkerrmat = std::nullopt;
8087 const AmgVector(5) & tpars = tmptrackpar->parameters();
8088 std::unique_ptr<const TrackParameters> trackpar(
8089 tmptrackpar->associatedSurface().createUniqueTrackParameters(tpars[0],
8094 std::move(trkerrmat))
8096 state->setTrackParameters(std::move(trackpar));
8097 FitQualityOnSurface fitQual{};
8099 if (errorok && trajectory.nDOF() > 0) {
8100 fitQual =
m_updator->fullStateFitQuality(
8101 *state->trackParameters(),
8102 measurement->localParameters(),
8103 measurement->localCovariance()
8106 fitQual = FitQualityOnSurface(0, state->numberOfMeasuredParameters());
8109 state->setFitQuality(fitQual);
8111 prevstate = state.get();
◆ calculateTrackParameters()
Definition at line 7598 of file GlobalChi2Fitter.cxx.
7606 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7607 int nstatesupstream = trajectory.numberOfUpstreamStates();
7608 const TrackParameters *prevtrackpar = trajectory.referenceParameters();
7609 std::unique_ptr<const TrackParameters> tmptrackpar;
7611 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7615 DistanceSolution distsol = surf1.straightLineDistanceEstimate(
7616 prevtrackpar->position(), prevtrackpar->momentum().unit()
7623 distsol.numberOfSolutions() > 0 &&
7624 prevtrackpar != trajectory.referenceParameters()
7634 trajectory.m_fieldprop,
7641 (
rv.m_parameters !=
nullptr) &&
7642 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7648 if (
rv.m_parameters ==
nullptr) {
7649 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7650 " pos: " << prevtrackpar->position() <<
" destination surface: " << surf1);
7654 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7658 if (
rv.m_jacobian != std::nullopt) {
7660 states[hitno]->materialEffects() !=
nullptr &&
7661 states[hitno]->materialEffects()->deltaE() != 0 &&
7662 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7663 !trajectory.m_straightline
7665 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7666 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7667 double mass = trajectory.mass();
7668 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7669 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7672 states[hitno]->setJacobian(*
rv.m_jacobian);
7673 }
else if (calcderiv) {
7678 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7680 if (meff !=
nullptr && hitno != 0) {
7681 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7682 surf, *meff, *
states[hitno]->trackParameters(), trajectory.mass(), -1
7685 if (std::holds_alternative<FitterStatusCode>(
r)) {
7686 return std::get<FitterStatusCode>(
r);
7689 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7690 prevtrackpar = tmptrackpar.get();
7692 prevtrackpar = currenttrackpar;
7696 prevtrackpar = trajectory.referenceParameters();
7698 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7701 DistanceSolution distsol = surf.straightLineDistanceEstimate(prevtrackpar->position(), prevtrackpar->momentum().unit());
7705 if (
distance < 0 && distsol.numberOfSolutions() > 0 && prevtrackpar != trajectory.referenceParameters()) {
7714 trajectory.m_fieldprop,
7720 (
rv.m_parameters !=
nullptr) &&
7722 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7727 if (
rv.m_parameters ==
nullptr) {
7728 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7729 " pos: " << prevtrackpar->
7730 position() <<
" destination surface: " << surf);
7734 if (
rv.m_jacobian != std::nullopt) {
7736 states[hitno]->materialEffects() !=
nullptr &&
7737 states[hitno]->materialEffects()->deltaE() != 0 &&
7738 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7739 !trajectory.m_straightline
7741 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7742 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7743 double mass = trajectory.mass();
7744 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7747 newp = std::sqrt(newp);
7750 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7753 states[hitno]->setJacobian(*
rv.m_jacobian);
7754 }
else if (calcderiv) {
7759 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7761 if (meff !=
nullptr) {
7762 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7763 surf, *meff, *
rv.m_parameters, trajectory.mass(), +1
7766 if (std::holds_alternative<FitterStatusCode>(
r)) {
7767 return std::get<FitterStatusCode>(
r);
7770 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7773 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7774 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 7572 of file GlobalChi2Fitter.cxx.
7581 PropagationResult
rv;
7584 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7587 if (
rv.m_parameters ==
nullptr) {
7588 propdir = invertPropdir(propdir);
7591 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 7533 of file GlobalChi2Fitter.cxx.
7542 std::unique_ptr<const TrackParameters>
rv;
7543 std::optional<TransportJacobian> jac{};
7554 if (
rv !=
nullptr && calcderiv) {
7559 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7565 return PropagationResult {
7568 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 5503 of file GlobalChi2Fitter.cxx.
5509 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5513 int nscatupstream = trajectory.numberOfUpstreamScatterers();
5514 int nbremupstream = trajectory.numberOfUpstreamBrems();
5515 int nscat = trajectory.numberOfScatterers();
5516 int nbrem = trajectory.numberOfBrems();
5517 int nperparams = trajectory.numberOfPerigeeParameters();
5519 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5522 int nmeas = (
int) weightderiv.rows();
5524 ParamDefsAccessor paraccessor;
5526 for (std::unique_ptr<GXFTrackState> & state :
states) {
5529 ((state->materialEffects() ==
nullptr) || state->materialEffects()->sigmaDeltaE() <= 0)
5535 const MeasurementBase *measbase = state->measurement();
5536 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5537 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5540 double sinstereo = 0;
5543 sinstereo = state->sinStereo();
5546 double cosstereo = (sinstereo == 0) ? 1. : std::sqrt(1 - sinstereo * sinstereo);
5548 for (
int i = 0;
i < 5;
i++) {
5550 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5556 if (trajectory.numberOfPerigeeParameters() > 0) {
5557 int cols = trajectory.m_straightline ? 4 : 5;
5560 weightderiv.row(measno).head(
cols) =
5561 (derivatives.row(0).head(
cols) * cosstereo +
5562 sinstereo * derivatives.row(1).head(
cols)) /
5565 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5569 for (
int j = scatmin; j < scatmax; j++) {
5570 int index = nperparams + ((trajectory.prefit() != 1) ? 2 * j : j);
5571 double thisderiv = 0;
5575 if (
i == 0 && sinstereo != 0) {
5576 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5581 weightderiv(measno,
index) = thisderiv /
error[measno];
5583 if (trajectory.prefit() != 1) {
5586 if (
i == 0 && sinstereo != 0) {
5587 thisderiv =
sign * (derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index));
5592 weightderiv(measno,
index) = thisderiv /
error[measno];
5596 for (
int j = bremmin; j < bremmax; j++) {
5597 double thisderiv = 0;
5598 int index = j + nperparams + 2 * nscat;
5599 if (
i == 0 && sinstereo != 0) {
5600 thisderiv = derivatives(0,
index) * cosstereo + sinstereo * derivatives(1,
index);
5602 thisderiv = derivatives(
i,
index);
5604 weightderiv(measno,
index) = thisderiv /
error[measno];
5610 double *
errors = state->measurementErrors();
5611 for (
int i = 0;
i < 5;
i++) {
5618 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5623 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5625 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5626 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5628 double mass = .001 * trajectory.mass();
5630 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5634 auto multiplier = [] (
double mass,
double qOverP){
5637 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5638 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5640 if (trajectory.numberOfPerigeeParameters() > 0) {
5641 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5644 const auto bremNoBase = nperparams + 2 * nscat;
5645 if (bremno < nbremupstream) {
5646 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5647 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5648 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5651 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5652 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5653 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
◆ fillResiduals()
Definition at line 5193 of file GlobalChi2Fitter.cxx.
5205 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5210 int nbrem = trajectory.numberOfBrems();
5211 int nperpars = trajectory.numberOfPerigeeParameters();
5212 int nfitpars = trajectory.numberOfFitParameters();
5215 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5216 int nidhits = trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits();
5217 int nsihits = trajectory.numberOfSiliconHits();
5218 int ntrthits = trajectory.numberOfTRTHits();
5219 int nhits = trajectory.numberOfHits();
5220 int nmeas = (
int)
res.size();
5222 ParamDefsAccessor paraccessor;
5223 bool scatwasupdated =
false;
5225 GXFTrackState *state_maxbrempull =
nullptr;
5226 int bremno_maxbrempull = 0;
5227 double maxbrempull = 0;
5229 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5230 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5233 const MeasurementBase *measbase = state->measurement();
5240 trajectory.nDOF() != 0 &&
5241 std::abs((trajectory.prevchi2() - trajectory.chi2()) / trajectory.nDOF()) < 15 &&
5242 !state->associatedSurface().isFree() &&
5243 nidhits < trajectory.numberOfHits() &&
5244 (nperpars == 0 || nidhits > 0) &&
5245 (!state->isRecalibrated())
5250 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5253 currenttrackpar->associatedSurface()
5256 state->setMeasurement(std::move(newpseudo));
5257 measbase = state->measurement();
5260 double *
errors = state->measurementErrors();
5264 for (
int i = 0;
i < 5;
i++) {
5266 !measbase->localParameters().contains(paraccessor.pardef[
i]) ||
5277 res[measno] = residuals[
i];
5279 if (
i == 2 && std::abs(std::abs(
res[measno]) - 2 *
M_PI) < std::abs(
res[measno])) {
5280 if (
res[measno] < 0) {
5289 double *
errors = state->measurementErrors();
5290 for (
int i = 0;
i < 5;
i++) {
5300 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5302 double deltaphi = state->materialEffects()->deltaPhi();
5303 double measdeltaphi = state->materialEffects()->measuredDeltaPhi();
5304 double sigmadeltaphi = state->materialEffects()->sigmaDeltaPhi();
5305 double deltatheta = state->materialEffects()->deltaTheta();
5306 double sigmadeltatheta = state->materialEffects()->sigmaDeltaTheta();
5308 if (trajectory.prefit() != 1) {
5309 b[nperpars + 2 * scatno] -= (deltaphi - measdeltaphi) / (sigmadeltaphi * sigmadeltaphi);
5310 b[nperpars + 2 * scatno + 1] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5312 b[nperpars + scatno] -= deltatheta / (sigmadeltatheta * sigmadeltatheta);
5316 deltaphi * deltaphi / (sigmadeltaphi * sigmadeltaphi) +
5317 deltatheta * deltatheta / (sigmadeltatheta * sigmadeltatheta)
5323 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5324 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5326 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5327 const double pbrem = 1. / std::abs(qoverpbrem);
5328 const double p = 1. / std::abs(qoverp);
5329 const double mass = .001 * trajectory.mass();
5331 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5333 res[nmeas - nbrem + bremno] = .001 * averagenergyloss -
energy + bremEnergy;
5335 double sigde = state->materialEffects()->sigmaDeltaE();
5336 double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5337 double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5339 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5341 if (state->materialEffects()->isKink()) {
5342 maxbrempull = -999999999;
5343 state_maxbrempull =
nullptr;
5347 cache.m_asymeloss &&
5349 (trajectory.prefit() == 0) &&
5351 sigde != sigdepos &&
5354 double elosspull =
res[nmeas - nbrem + bremno] / (.001 * sigde);
5356 if (trajectory.mass() > 100) {
5357 if (elosspull < -1) {
5358 state->materialEffects()->setSigmaDeltaE(sigdepos);
5359 }
else if (elosspull > 1) {
5360 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5363 error[nmeas - nbrem + bremno] = .001 * state->materialEffects()->sigmaDeltaE();
5364 }
else if ((trajectory.numberOfTRTHits() == 0) ||
it >= 3) {
5366 !state->materialEffects()->isKink() && (
5367 (elosspull < -.2 &&
m_fixbrem == -1 && elosspull < maxbrempull) ||
5371 bremno_maxbrempull = bremno;
5372 state_maxbrempull = state.get();
5373 maxbrempull = elosspull;
5382 (trajectory.prefit() == 0) &&
5383 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5384 state->materialEffects()->isMeasuredEloss() &&
5385 res[nmeas - nbrem + bremno] / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5388 trajectory.prefit() != 0 ? trajectory.referenceParameters()
5389 :
states[hitno - 2]->trackParameters();
5392 std::vector<MaterialEffectsOnTrack> calomeots =
5397 parforcalo->associatedSurface(),
5401 if (calomeots.size() == 3) {
5402 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5403 double newres = .001 * averagenergyloss -
energy + bremEnergy;
5404 double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5406 if (std::abs(newres / newerr) < std::abs(
res[nmeas - nbrem + bremno] /
error[nmeas - nbrem + bremno])) {
5407 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5409 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5410 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5411 res[nmeas - nbrem + bremno] = newres;
5412 error[nmeas - nbrem + bremno] = newerr;
5416 state->materialEffects()->setMeasuredEloss(
false);
5424 for (; measno < nmeas; measno++) {
5425 if (
error[measno] == 0) {
5432 if (!doderiv && (scatwasupdated)) {
5436 double oldchi2 = trajectory.chi2();
5437 trajectory.setPrevChi2(oldchi2);
5438 trajectory.setChi2(
chi2);
5440 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5441 double newredchi2 = (trajectory.nDOF() > 0) ?
chi2 / trajectory.nDOF() : 0;
5444 trajectory.prefit() > 0 && (
5445 (newredchi2 < 2 &&
it != 0) ||
5446 (newredchi2 < oldredchi2 + .1 && std::abs(newredchi2 - oldredchi2) < 1 &&
it != 1)
5449 trajectory.setConverged(
true);
5452 double maxdiff = (nsihits != 0 && nsihits + ntrthits == nhits &&
chi2 < oldchi2) ? 200 : 1.;
5454 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5456 if (miniter < cache.m_miniter) {
5457 miniter = cache.m_miniter;
5461 trajectory.setConverged(
true);
5464 if ((state_maxbrempull !=
nullptr) && trajectory.converged()) {
5465 state_maxbrempull->materialEffects()->setSigmaDeltaE(
5466 10 * state_maxbrempull->materialEffects()->sigmaDeltaEPos()
5469 state_maxbrempull->materialEffects()->setKink(
true);
5470 trajectory.setConverged(
false);
5472 double olderror =
error[nmeas - nbrem + bremno_maxbrempull];
5473 double newerror = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5474 error[nmeas - nbrem + bremno_maxbrempull] = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5476 if (
a.cols() != nfitpars) {
5480 for (
int i = 0;
i < nfitpars;
i++) {
5481 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5485 for (
int j =
i; j < nfitpars; j++) {
5489 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5490 weightderiv(nmeas - nbrem + bremno_maxbrempull, j) *
5491 (1 - olderror * olderror / (newerror * newerror))
5495 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= olderror / newerror;
5498 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 7482 of file GlobalChi2Fitter.cxx.
7492 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7493 ctx,
src, dst.associatedSurface(), propdir,
false
7505 &
rv.front()->associatedSurface() == &dst.associatedSurface() ||
7506 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7507 trackParametersClose(*
rv.front(),
src, 0.001) ||
7508 trackParametersClose(*
rv.front(), *dst.trackParameters(), 0.001)
7511 rv.front().reset(
nullptr);
7521 &
rv.back()->associatedSurface() == &dst.associatedSurface() ||
7522 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7523 trackParametersClose(*
rv.back(),
src, 0.001) ||
7524 trackParametersClose(*
rv.back(), *dst.trackParameters(), 0.001)
7527 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 7033 of file GlobalChi2Fitter.cxx.
7046 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7052 if (
tp ==
nullptr) {
7061 const TrkDetElementBase * de =
tp->associatedSurface().associatedDetectorElement();
7063 if (de ==
nullptr) {
7074 if (id_set.find(
id) != id_set.end()) {
7123 if (sct_set.find(
os) != sct_set.end()) {
7124 ++
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 7218 of file GlobalChi2Fitter.cxx.
7232 constexpr
uint min_meas = 3;
7233 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7237 bool seen_meas =
false;
7239 std::set<Identifier> id_set;
7240 std::set<Identifier> sct_set;
7246 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7269 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7276 if (seen_meas && dist >= 2.5) {
7283 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7284 std::vector<std::unique_ptr<TrackParameters>>
states;
7292 if (hc.has_value()) {
7313 GXFTrackState & last =
states.back();
7324 last.trackParameters() !=
nullptr &&
7331 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7333 *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 7143 of file GlobalChi2Fitter.cxx.
7151 GXFTrackState * lastmeas =
nullptr;
7153 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7165 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7166 rv.reserve(trajectory.trackStates().size());
7173 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7187 rv.emplace_back(*
s);
7194 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7196 if (de !=
nullptr) {
7209 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 8349 of file GlobalChi2Fitter.cxx.
8359 if (cond_obj ==
nullptr) {
8364 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 8307 of file GlobalChi2Fitter.cxx.
8308 const auto *pDataVector = intrk1.measurementsOnTrack();
8309 auto nmeas1 = pDataVector->size();
8310 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8318 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8320 if (lastMeasIsCompetingRIO){
8322 testrot = &testcrot->rioOnTrack(0);
8326 if (testrot ==
nullptr) {
8327 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8330 if(penultimateIsRIO){
8331 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8333 if (penultimateIsCompetingRIO){
8335 testrot = &testcrot->rioOnTrack(0);
8342 (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 7351 of file GlobalChi2Fitter.cxx.
7358 auto trajectory = std::make_unique<Trk::TrackStates>();
7364 GXFTrajectory tmptrajectory(oldtrajectory);
7366 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7368 if (perigee_ts ==
nullptr) {
7372 tmptrajectory.addBasicState(std::move(perigee_ts), cache.m_acceleration ? 0 : tmptrajectory.numberOfUpstreamStates());
7374 trajectory->reserve(tmptrajectory.trackStates().size());
7375 for (
auto & hit : tmptrajectory.trackStates()) {
7380 hit->resetTrackCovariance();
7386 hit->materialEffects()))
7392 auto trackState = hit->trackStateOnSurface();
7393 hit->resetTrackCovariance();
7394 trajectory->emplace_back(trackState.release());
7397 auto qual = std::make_unique<FitQuality>(tmptrajectory.chi2(), tmptrajectory.nDOF());
7408 if (matEffects ==
electron && tmptrajectory.hasKink()) {
7413 if (tmptrajectory.m_straightline) {
7417 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7427 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7436 std::optional<TrackHoleCount> hole_count;
7461 if (hole_count.has_value()) {
7474 rv->setTrackSummary(std::move(
ts));
◆ makeTrackFillDerivativeMatrix()
void Trk::GlobalChi2Fitter::makeTrackFillDerivativeMatrix |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
oldtrajectory |
|
) |
| |
|
staticprivate |
Definition at line 6815 of file GlobalChi2Fitter.cxx.
6819 Amg::MatrixX & derivs = oldtrajectory.weightedResidualDerivatives();
6823 for (
auto & hit : oldtrajectory.trackStates()) {
6824 if (
const auto *pMeas{hit->measurement()};
6830 nrealmeas += hit->numberOfMeasuredParameters();
6833 cache.m_derivmat.resize(nrealmeas, oldtrajectory.numberOfFitParameters());
6834 cache.m_derivmat.setZero();
6837 int nperpars = oldtrajectory.numberOfPerigeeParameters();
6838 int nscat = oldtrajectory.numberOfScatterers();
6839 for (
auto & hit : oldtrajectory.trackStates()) {
6840 if (
const auto *pMeas{hit->measurement()};
6846 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
6847 for (
int j = 0; j < oldtrajectory.numberOfFitParameters(); j++) {
6848 cache.m_derivmat(
i, j) = derivs(measindex2, j) *
errors[measindex2];
6849 if ((j == 4 && !oldtrajectory.m_straightline) || j >= nperpars + 2 * nscat) {
6850 cache.m_derivmat(
i, j) *= 1000;
6857 measindex += hit->numberOfMeasuredParameters();
6858 }
else if (hit->materialEffects() ==
nullptr) {
6859 measindex2 += hit->numberOfMeasuredParameters();
◆ makeTrackFindPerigee()
◆ makeTrackFindPerigeeParameters()
Definition at line 6864 of file GlobalChi2Fitter.cxx.
6870 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
6871 std::tie(firstmeasstate, lastmeasstate) = oldtrajectory.findFirstLastMeasurement();
6872 std::unique_ptr<const TrackParameters> per(
nullptr);
6875 std::unique_ptr<const TrackParameters> prevpar(
6876 firstmeasstate->trackParameters() !=
nullptr ?
6877 firstmeasstate->trackParameters()->clone() :
6880 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.upstreamMaterialLayers();
6883 for (
int i = (
int)upstreamlayers.size() - 1;
i >= 0;
i--) {
6884 if (prevpar ==
nullptr) {
6891 if (
layer ==
nullptr) {
6892 layer = upstreamlayers[
i].second;
6895 DistanceSolution distsol =
layer->surfaceRepresentation().straightLineDistanceEstimate(
6896 prevpar->position(), prevpar->momentum().unit()
6900 if (distsol.numberOfSolutions() == 2) {
6905 if (distsol.first() * distsol.second() < 0 && !
first) {
6914 std::unique_ptr<const TrackParameters> layerpar(
6918 layer->surfaceRepresentation(),
6921 oldtrajectory.m_fieldprop,
6926 if (layerpar ==
nullptr) {
6930 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
6934 prevpar = std::move(layerpar);
6938 const Layer *startlayer = firstmeasstate->trackParameters()->associatedSurface().associatedLayer();
6940 if ((startlayer !=
nullptr) && (startlayer->layerMaterialProperties() !=
nullptr)) {
6941 double startfactor = startlayer->layerMaterialProperties()->alongPostFactor();
6942 const Surface & discsurf = startlayer->surfaceRepresentation();
6945 startfactor = startlayer->layerMaterialProperties()->oppositePostFactor();
6947 if (startfactor > 0.5) {
6948 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6949 firstmeasstate->trackParameters(), *startlayer,
oppositeMomentum, matEffects
6952 if (updatedpar !=
nullptr) {
6953 firstmeasstate->setTrackParameters(std::move(updatedpar));
6962 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
6964 if ((endlayer !=
nullptr) && (endlayer->layerMaterialProperties() !=
nullptr)) {
6965 double endfactor = endlayer->layerMaterialProperties()->alongPreFactor();
6966 const Surface & discsurf = endlayer->surfaceRepresentation();
6969 endfactor = endlayer->layerMaterialProperties()->oppositePreFactor();
6972 if (endfactor > 0.5) {
6973 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
6974 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
6977 if (updatedpar !=
nullptr) {
6978 lastmeasstate->setTrackParameters(std::move(updatedpar));
6983 if (prevpar !=
nullptr) {
6990 oldtrajectory.m_fieldprop,
6995 if (per ==
nullptr) {
6996 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7001 }
else if (cache.m_acceleration && (firstmeasstate->trackParameters() !=
nullptr)) {
7003 *firstmeasstate->trackParameters(),
7009 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) {
5099 delete finaltrajectory;
5108 if (!cache.m_acceleration && (finaltrajectory->prefit() == 0)) {
5109 if (nperpars == 5) {
5110 for (
int i = 0;
i <
a.cols();
i++) {
5111 a_inv(4,
i) *= .001;
5112 a_inv(
i, 4) *= .001;
5116 int scatterPos = nperpars + 2 * nscat;
5117 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
5118 for (
int i = 0;
i <
a.cols();
i++) {
5119 a_inv(scatterPos,
i) *= .001;
5120 a_inv(
i, scatterPos) *= .001;
5126 int nperparams = finaltrajectory->numberOfPerigeeParameters();
5127 for (
int i = 0;
i < nperparams;
i++) {
5128 for (
int j = 0; j < nperparams; j++) {
5129 (errmat) (j,
i) = a_inv(j,
i);
5133 if (trajectory.m_straightline) {
5134 (errmat) (4, 4) = 1
e-20;
5137 const AmgVector(5) & perpars = finaltrajectory->referenceParameters()->parameters();
5138 std::unique_ptr<const TrackParameters> measper(
5139 finaltrajectory->referenceParameters()->associatedSurface().createUniqueTrackParameters(
5140 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5144 finaltrajectory->setReferenceParameters(std::move(measper));
5146 cache.m_fullcovmat = a_inv;
5150 std::unique_ptr<Track>
track =
nullptr;
5152 if (finaltrajectory->prefit() > 0) {
5153 if (finaltrajectory != &trajectory) {
5154 delete finaltrajectory;
5159 if (finaltrajectory->numberOfOutliers() <=
m_maxoutliers || !runOutlier) {
5166 double cut = (finaltrajectory->numberOfSiliconHits() ==
5167 finaltrajectory->numberOfHits())
5173 (
track !=
nullptr) && (
5174 track->fitQuality()->numberDoF() != 0 &&
5175 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5178 track.reset(
nullptr);
5182 if (
track ==
nullptr) {
5186 if (finaltrajectory != &trajectory) {
5187 delete finaltrajectory;
5190 return track.release();
◆ myfit_helper()
◆ numericalDerivatives()
Definition at line 8116 of file GlobalChi2Fitter.cxx.
8123 ParamDefsAccessor paraccessor;
8131 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8134 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8141 const Surface & previousSurface = tmpprevpar->associatedSurface();
8145 for (
int i = 0;
i < 5;
i++) {
8148 if (thisdiscsurf &&
i == 1) {
8152 vecpluseps[paraccessor.pardef[
i]] += eps[
i];
8153 vecminuseps[paraccessor.pardef[
i]] -= eps[
i];
8154 if (thiscylsurf &&
i == 0) {
8155 if (vecpluseps[0] / previousSurface.bounds().r() >
M_PI) {
8156 vecpluseps[0] -= 2 *
M_PI * previousSurface.bounds().r();
8158 if (vecminuseps[0] / previousSurface.bounds().r() < -
M_PI) {
8159 vecminuseps[0] += 2 *
M_PI * previousSurface.bounds().r();
8162 if (thisdiscsurf &&
i == 1) {
8163 if (vecpluseps[
i] >
M_PI) {
8164 vecpluseps[
i] -= 2 *
M_PI;
8166 if (vecminuseps[
i] < -
M_PI) {
8167 vecminuseps[
i] += 2 *
M_PI;
8173 std::unique_ptr<const TrackParameters> parpluseps(
8174 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8183 std::unique_ptr<const TrackParameters> parminuseps(
8184 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8194 std::unique_ptr<const TrackParameters> newparpluseps(
8205 std::unique_ptr<const TrackParameters> newparminuseps(
8220 if (newparpluseps ==
nullptr) {
8232 if (newparminuseps ==
nullptr) {
8244 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8248 for (
int j = 0; j < 5; j++) {
8249 double diff = newparpluseps->parameters()[paraccessor.pardef[j]] -
8250 newparminuseps->parameters()[paraccessor.pardef[j]];
8251 if (cylsurf && j == 0) {
8252 double length = 2 *
M_PI * surf.bounds().r();
8261 if (discsurf && j == 1) {
8262 if (std::abs(std::abs(
diff) - 2 *
M_PI) < std::abs(
diff)) {
8271 (*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 5661 of file GlobalChi2Fitter.cxx.
5672 int nfitpars = trajectory.numberOfFitParameters();
5673 int nperpars = trajectory.numberOfPerigeeParameters();
5674 int scatpars = 2 * trajectory.numberOfScatterers();
5675 int nupstreamstates = trajectory.numberOfUpstreamStates();
5676 int nbrem = trajectory.numberOfBrems();
5677 double oldchi2 = trajectory.chi2();
5678 double oldredchi2 = (trajectory.nDOF() > 0) ? oldchi2 / trajectory.nDOF() : 0;
5679 int nsihits = trajectory.numberOfSiliconHits();
5680 int ntrthits = trajectory.numberOfTRTHits();
5681 int nhits = trajectory.numberOfHits();
5683 if (cache.m_phiweight.empty()) {
5684 cache.m_phiweight.assign(trajectory.trackStates().size(), 1);
5697 double newredchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
5699 ATH_MSG_DEBUG(
"old chi2: " << oldchi2 <<
"/" << trajectory.nDOF() <<
5700 "=" << oldredchi2 <<
" new chi2: " << trajectory.chi2() <<
"/" <<
5701 trajectory.nDOF() <<
"=" << newredchi2);
5703 if (trajectory.prefit() > 0 && trajectory.converged()) {
5709 std::vector < std::pair < double, double >>&scatsigmas = trajectory.scatteringSigmas();
5711 int nmeas = (
int)
res.size();
5713 const Amg::MatrixX & weight_deriv = trajectory.weightedResidualDerivatives();
5720 if (cache.m_firstmeasurement.empty()) {
5721 cache.m_firstmeasurement.resize(nfitpars);
5722 cache.m_lastmeasurement.resize(nfitpars);
5723 for (
int i = 0;
i < nperpars;
i++) {
5724 cache.m_firstmeasurement[
i] = 0;
5725 cache.m_lastmeasurement[
i] = nmeas - nbrem;
5730 for (
int i = 0;
i < (
int) trajectory.trackStates().size();
i++) {
5731 std::unique_ptr<GXFTrackState> & state = trajectory.trackStates()[
i];
5732 GXFMaterialEffects *meff = state->materialEffects();
5733 if (meff ==
nullptr) {
5734 measno += state->numberOfMeasuredParameters();
5736 if (meff !=
nullptr) {
5737 if (meff->sigmaDeltaTheta() != 0
5738 && ((trajectory.prefit() == 0) || meff->deltaE() == 0)) {
5739 int scatterPos = nperpars + 2 * scatno;
5740 if (
i < nupstreamstates) {
5741 cache.m_lastmeasurement[scatterPos] =
5742 cache.m_lastmeasurement[scatterPos + 1] = measno;
5743 cache.m_firstmeasurement[scatterPos] =
5744 cache.m_firstmeasurement[scatterPos + 1] = 0;
5746 cache.m_lastmeasurement[scatterPos] =
5747 cache.m_lastmeasurement[scatterPos + 1] = nmeas - nbrem;
5748 cache.m_firstmeasurement[scatterPos] =
5749 cache.m_firstmeasurement[scatterPos + 1] = measno;
5753 if (meff->sigmaDeltaE() > 0) {
5754 if (
i < nupstreamstates) {
5755 cache.m_firstmeasurement[nperpars + scatpars + bremno] = 0;
5756 cache.m_lastmeasurement[nperpars + scatpars + bremno] = measno;
5758 cache.m_firstmeasurement[nperpars + scatpars + bremno] = measno;
5759 cache.m_lastmeasurement[nperpars + scatpars + bremno] =
5769 if (
a.cols() != nfitpars) {
5773 for (
int k = 0;
k < nfitpars;
k++) {
5775 int maxmeas = nmeas - nbrem;
5776 maxmeas = cache.m_lastmeasurement[
k];
5777 minmeas = cache.m_firstmeasurement[
k];
5779 for (measno = minmeas; measno < maxmeas; measno++) {
5781 res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5785 if (
k == 4 ||
k >= nperpars + scatpars) {
5786 for (measno = nmeas - nbrem; measno < nmeas; measno++) {
5787 b[
k] +=
res[measno] * (1. /
error[measno]) * weight_deriv(measno,
k);
5792 for (
int l =
k;
l < nfitpars;
l++) {
5794 std::min(cache.m_lastmeasurement[
k], cache.m_lastmeasurement[
l]);
5797 cache.m_firstmeasurement[
l]);
5799 for (measno = minmeas; measno < maxmeas; measno++) {
5800 tmp += weight_deriv(measno,
k) * weight_deriv(measno,
l);
5802 a.fillSymmetric(
l,
k,
tmp);
5810 for (
int k = nperpars;
k < nperpars + scatpars;
k += 2) {
5811 a(
k,
k) += 1. / (scatsigmas[scatno].first * scatsigmas[scatno].first);
5812 a(
k + 1,
k + 1) += 1. / (scatsigmas[scatno].second * scatsigmas[scatno].second);
5816 for (
int measno = nmeas - nbrem; measno < nmeas; measno++) {
5817 for (
int k = 4;
k < nfitpars;
k++) {
5819 k = nperpars + scatpars;
5822 for (
int l =
k;
l < nfitpars;
l++) {
5824 l = nperpars + scatpars;
5826 double tmp =
a(
l,
k) + weight_deriv(measno,
k) * weight_deriv(measno,
l);
5827 a.fillSymmetric(
l,
k,
tmp);
5833 unsigned int scatno = 0;
5834 bool weightchanged =
false;
5836 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5837 GXFMaterialEffects *meff = thisstate->materialEffects();
5839 if (meff !=
nullptr) {
5840 const PlaneSurface *plsurf =
nullptr;
5843 plsurf =
static_cast < const PlaneSurface *
>(&thisstate->associatedSurface());
5844 if (meff->deltaE() == 0 || ((trajectory.prefit() == 0) && (plsurf !=
nullptr))) {
5845 weightchanged =
true;
5847 if (
a.cols() != nfitpars) {
5851 int scatNoIndex = 2 * scatno + nperpars;
5853 if (trajectory.prefit() == 0) {
5854 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5855 if (scatno >= cache.m_phiweight.size()) {
5857 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5858 throw std::range_error(
message.str());
5862 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5866 cache.m_phiweight[scatno] = 1.00000001;
5867 }
else if (
it == 1) {
5868 cache.m_phiweight[scatno] = 1.0000001;
5869 }
else if (
it <= 3) {
5870 cache.m_phiweight[scatno] = 1.0001;
5871 }
else if (
it <= 6) {
5872 cache.m_phiweight[scatno] = 1.01;
5874 cache.m_phiweight[scatno] = 1.1;
5877 a(scatNoIndex, scatNoIndex) *= cache.m_phiweight[scatno];
5881 else if (trajectory.prefit() >= 2) {
5882 if (newredchi2 > oldredchi2 - 1 && newredchi2 < oldredchi2) {
5883 a(scatNoIndex, scatNoIndex) *= 1.0001;
5884 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5885 }
else if (newredchi2 > oldredchi2 - 25 && newredchi2 < oldredchi2) {
5886 a(scatNoIndex, scatNoIndex) *= 1.001;
5887 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.0001;
5889 a(scatNoIndex, scatNoIndex) *= 1.1;
5890 a(scatNoIndex + 1, scatNoIndex + 1) *= 1.001;
5896 thisstate->materialEffects()->sigmaDeltaPhi() != 0 &&
5897 ((trajectory.prefit() == 0) || thisstate->materialEffects()->deltaE() == 0)
5905 (trajectory.prefit() == 2) &&
5907 trajectory.numberOfBrems() > 0 &&
5908 (newredchi2 < oldredchi2 - 25 || newredchi2 > oldredchi2)
5913 if (doderiv || weightchanged) {
5917 if (trajectory.converged()) {
5918 if ((trajectory.prefit() == 0) && nsihits + ntrthits != nhits) {
5919 unsigned int scatno = 0;
5921 if (
a.cols() != nfitpars) {
5925 for (std::unique_ptr<GXFTrackState> & thisstate : trajectory.trackStates()) {
5926 if ((thisstate->materialEffects() !=
nullptr) && thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5927 if (scatno >= cache.m_phiweight.size()) {
5929 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5930 throw std::range_error(
message.str());
5933 const PlaneSurface *plsurf =
nullptr;
5936 plsurf =
static_cast<const PlaneSurface *
>(&thisstate->associatedSurface());
5938 if (thisstate->materialEffects()->deltaE() == 0 || (plsurf !=
nullptr)) {
5939 int scatNoIndex = 2 * scatno + nperpars;
5940 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5941 cache.m_phiweight[scatno] = 1;
5944 if (thisstate->materialEffects()->sigmaDeltaPhi() != 0) {
5957 (newredchi2 < 2 || (newredchi2 < oldredchi2 && newredchi2 > oldredchi2 - .5)) &&
5958 (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 6321 of file GlobalChi2Fitter.cxx.
6330 bool trackok =
false;
6331 GXFTrajectory *oldtrajectory = &trajectory;
6332 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6333 GXFTrajectory *newtrajectory =
nullptr;
6334 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6340 while (!trackok && oldtrajectory->nDOF() > 0) {
6342 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->trackStates();
6345 Amg::MatrixX & weightderiv = oldtrajectory->weightedResidualDerivatives();
6346 int nfitpars = oldtrajectory->numberOfFitParameters();
6347 int nhits = oldtrajectory->numberOfHits();
6348 int nsihits = oldtrajectory->numberOfSiliconHits();
6350 if (nhits != nsihits) {
6354 double maxsipull = -1;
6356 int hitno_maxsipull = -1;
6357 int measno_maxsipull = -1;
6358 int stateno_maxsipull = 0;
6359 GXFTrackState *state_maxsipull =
nullptr;
6364 int noutl = oldtrajectory->numberOfOutliers();
6370 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6371 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6377 double *
errors = state->measurementErrors();
6379 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6380 double sinstereo = state->sinStereo();
6381 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6382 double weight1 = -1;
6384 if (hitcov(0, 0) > trackcov(0, 0)) {
6385 if (sinstereo == 0) {
6389 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6390 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6401 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6404 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6407 sipull1 =
std::max(sipull1, sipull2);
6409 if (sipull1 > maxsipull) {
6410 maxsipull = sipull1;
6411 measno_maxsipull = measno;
6412 state_maxsipull = state.get();
6413 stateno_maxsipull = stateno;
6414 hitno_maxsipull = hitno;
6425 measno += state->numberOfMeasuredParameters();
6429 double maxpull = maxsipull;
6431 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6432 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6441 oldtrajectory->chi2() / oldtrajectory->nDOF() > .25 *
m_chi2cut
6443 state_maxsipull = oldtrajectory->trackStates()[stateno_maxsipull].get();
6444 const PrepRawData *prd{};
6446 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6447 prd = rot->prepRawData();
6449 std::unique_ptr < const RIO_OnTrack > broadrot;
6450 double *olderror = state_maxsipull->measurementErrors();
6452 const TrackParameters *trackpar_maxsipull = state_maxsipull->trackParameters();
6454 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6455 std::unique_ptr<const TrackParameters> trackparForCorrect(
6456 trackpar_maxsipull->associatedSurface().createUniqueTrackParameters(
6462 state_maxsipull->hasTrackCovariance()
6464 state_maxsipull->trackCovariance())
6468 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6469 double newpull = -1;
6470 double newpull1 = -1;
6471 double newpull2 = -1;
6472 double newres1 = -1;
6473 double newres2 = -1;
6474 double newsinstereo = 0;
6478 !state_maxsipull->isRecalibrated() &&
6480 oldtrajectory->chi2() / trajectory.nDOF() > .3 *
m_chi2cut &&
6488 newerror[0] = std::sqrt(
covmat(0, 0));
6490 if (state_maxsipull->sinStereo() != 0) {
6508 newerror[0] = std::sqrt(
v0);
6511 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6513 if (cosstereo != 1.) {
6515 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6516 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6519 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6522 if (newerror[0] == 0.0) {
6523 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6527 newpull1 = std::abs(newres1 / newerror[0]);
6531 newerror[1] = std::sqrt(
covmat(1, 1));
6532 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6533 newpull2 = std::abs(newres2 / newerror[1]);
6536 newpull =
std::max(newpull1, newpull2);
6542 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6544 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6545 throw std::runtime_error(
6546 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6551 newtrajectory = oldtrajectory;
6553 if (
a.cols() != nfitpars) {
6557 double oldres1 =
res[measno_maxsipull];
6558 res[measno_maxsipull] = newres1;
6559 err[measno_maxsipull] = newerror[0];
6561 for (
int i = 0;
i < nfitpars;
i++) {
6562 if (weightderiv(measno_maxsipull,
i) == 0) {
6566 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6568 for (
int j =
i; j < nfitpars; j++) {
6572 weightderiv(measno_maxsipull,
i) *
6573 weightderiv(measno_maxsipull, j) *
6574 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6578 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6582 double oldres2 =
res[measno_maxsipull + 1];
6583 res[measno_maxsipull + 1] = newres2;
6584 err[measno_maxsipull + 1] = newerror[1];
6586 for (
int i = 0;
i < nfitpars;
i++) {
6587 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6591 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6593 for (
int j =
i; j < nfitpars; j++) {
6597 weightderiv(measno_maxsipull + 1,
i) *
6598 weightderiv(measno_maxsipull + 1, j) *
6599 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6604 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6609 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6610 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6611 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6612 olderror[1] <<
" newerror_1=" << newerror[1]
6615 state_maxsipull->setMeasurement(std::move(broadrot));
6616 state_maxsipull->setSinStereo(newsinstereo);
6617 state_maxsipull->setMeasurementErrors(newerror);
6621 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6622 (oldtrajectory->chi2() / oldtrajectory->nDOF() > .3 *
m_chi2cut || noutl > 1)
6626 (oldtrajectory->nDOF() > 1 || hittype_maxsipull ==
TrackState::SCT) &&
6631 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6632 measno_maxsipull <<
" pull=" << maxsipull
6639 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6640 newtrajectory = cleanup_newtrajectory.get();
6642 if (newa.cols() != nfitpars) {
6647 Amg::MatrixX & newweightderiv = newtrajectory->weightedResidualDerivatives();
6648 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6649 throw std::runtime_error(
6650 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6654 double oldres1 =
res[measno_maxsipull];
6655 newres[measno_maxsipull] = 0;
6657 for (
int i = 0;
i < nfitpars;
i++) {
6658 if (weightderiv(measno_maxsipull,
i) == 0) {
6662 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6664 for (
int j =
i; j < nfitpars; j++) {
6668 weightderiv(measno_maxsipull,
i) *
6669 weightderiv(measno_maxsipull, j)
6673 newweightderiv(measno_maxsipull,
i) = 0;
6677 double oldres2 =
res[measno_maxsipull + 1];
6678 newres[measno_maxsipull + 1] = 0;
6680 for (
int i = 0;
i < nfitpars;
i++) {
6681 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6685 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6687 for (
int j =
i; j < nfitpars; j++) {
6688 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6695 weightderiv(measno_maxsipull + 1,
i) *
6696 weightderiv(measno_maxsipull + 1, j)
6700 newweightderiv(measno_maxsipull + 1,
i) = 0;
6704 newtrajectory->setOutlier(stateno_maxsipull);
6710 newtrajectory->setConverged(
false);
6726 if (!newtrajectory->converged()) {
6728 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6735 if (!newtrajectory->converged()) {
6744 double oldchi2 = oldtrajectory->chi2() / oldtrajectory->nDOF();
6745 double newchi2 = (newtrajectory->nDOF() > 0) ? newtrajectory->chi2() / newtrajectory->nDOF() : 0;
6748 if (newtrajectory->nDOF() != oldtrajectory->nDOF() && maxsipull > cut2) {
6749 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6751 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6756 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6757 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6765 (void)cleanup_oldtrajectory.release();
6766 return oldtrajectory;
6768 if (oldtrajectory != newtrajectory) {
6769 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6770 oldtrajectory = newtrajectory;
6777 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6778 if (lltOfW.info() == Eigen::Success) {
6782 int ncols =
a.cols();
6783 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6784 fullcov = lltOfW.solve(weightInvAMG);
6802 oldtrajectory->nDOF() > 0 &&
6803 oldtrajectory->chi2() / oldtrajectory->nDOF() >
m_chi2cut &&
6811 (void)cleanup_oldtrajectory.release();
6812 return oldtrajectory;
◆ runTrackCleanerTRT()
Definition at line 6151 of file GlobalChi2Fitter.cxx.
6164 if (
it == 1 && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
6168 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6171 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6172 int nfitpars = trajectory.numberOfFitParameters();
6174 if (
a.cols() != nfitpars) {
6178 int nperpars = trajectory.numberOfPerigeeParameters();
6179 int nscats = trajectory.numberOfScatterers();
6182 bool outlierremoved =
false;
6183 bool hitrecalibrated =
false;
6185 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6186 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6194 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6198 trajectory.setOutlier(stateno);
6199 outlierremoved =
true;
6201 double *
errors = state->measurementErrors();
6202 double olderror =
errors[0];
6204 trajectory.updateTRTHitCount(stateno, olderror);
6206 for (
int i = 0;
i < nfitpars;
i++) {
6207 if (weightderiv(measno,
i) == 0) {
6211 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6213 for (
int j =
i; j < nfitpars; j++) {
6216 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6219 weightderiv(measno,
i) = 0;
6223 }
else if (trtrecal) {
6224 double *
errors = state->measurementErrors();
6225 double olderror =
errors[0];
6227 const auto *
const thisMeasurement{state->measurement()};
6233 if (oldrot->prepRawData() !=
nullptr) {
6234 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6236 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6238 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6239 double distance = std::abs(std::abs(trackradius) - dcradius);
6241 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6242 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6243 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6244 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6247 if (newrot !=
nullptr) {
6249 hitrecalibrated =
true;
6253 if ((measno < 0) or (measno >= (
int)
res.size())) {
6254 throw std::runtime_error(
6255 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6259 double oldres =
res[measno];
6260 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6262 state->setMeasurement(std::move(newrot));
6264 trajectory.updateTRTHitCount(stateno, olderror);
6266 for (
int i = 0;
i < nfitpars;
i++) {
6267 if (weightderiv(measno,
i) == 0) {
6271 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6273 for (
int j =
i; j < nfitpars; j++) {
6277 !cache.m_phiweight.empty() &&
6280 i < nperpars + 2 * nscats &&
6281 (
i - nperpars) % 2 == 0
6283 weight = cache.m_phiweight[(
i - nperpars) / 2];
6288 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6291 weightderiv(measno,
i) *= olderror / newerror;
6294 res[measno] = newres;
6295 err[measno] = newerror;
6304 measno += state->numberOfMeasuredParameters();
6308 if (trajectory.nDOF() < 0) {
6313 if (outlierremoved || hitrecalibrated) {
6315 trajectory.setConverged(
false);
6317 cache.m_miniter =
it + 2;
◆ setMinIterations()
void Trk::GlobalChi2Fitter::setMinIterations |
( |
int |
| ) |
|
|
privatevirtual |
Definition at line 8282 of file GlobalChi2Fitter.cxx.
8284 (
"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 7780 of file GlobalChi2Fitter.cxx.
7793 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7797 double newqoverp = 0;
7799 if (meff.sigmaDeltaE() <= 0) {
7804 double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.deltaE()) * std::sqrt(
mass *
mass + oldp * oldp) + meff.deltaE() * meff.deltaE();
7811 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
7817 return surf.createUniqueTrackParameters(
7818 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
◆ updateFitParameters()
Definition at line 5966 of file GlobalChi2Fitter.cxx.
5974 double d0 = refpar->parameters()[
Trk::d0];
5975 double z0 = refpar->parameters()[
Trk::z0];
5978 double qoverp = refpar->parameters()[
Trk::qOverP];
5979 int nscat = trajectory.numberOfScatterers();
5980 int nbrem = trajectory.numberOfBrems();
5981 int nperparams = trajectory.numberOfPerigeeParameters();
5983 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
5986 if (llt.info() == Eigen::Success) {
5992 if (trajectory.numberOfPerigeeParameters() > 0) {
5997 qoverp = (trajectory.m_straightline) ? 0 : .001 *
result[4] + qoverp;
6006 std::vector < std::pair < double, double >>&scatangles = trajectory.scatteringAngles();
6007 std::vector < double >&delta_ps = trajectory.brems();
6009 for (
int i = 0;
i < nscat;
i++) {
6010 scatangles[
i].first +=
result[2 *
i + nperparams];
6011 scatangles[
i].second +=
result[2 *
i + nperparams + 1];
6014 for (
int i = 0;
i < nbrem;
i++) {
6015 delta_ps[
i] +=
result[nperparams + 2 * nscat +
i];
6018 std::unique_ptr<const TrackParameters> newper(
6019 trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
6024 trajectory.setReferenceParameters(std::move(newper));
6025 trajectory.setScatteringAngles(scatangles);
6026 trajectory.setBrems(delta_ps);
◆ updatePixelROTs()
Update the Pixel ROT using the current trajectory/local track parameters.
Definition at line 6031 of file GlobalChi2Fitter.cxx.
6037 if ( trajectory.numberOfSiliconHits() == 0) {
6046 if (!splitProbContainer.isValid()) {
6050 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6053 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6054 int nfitpars = trajectory.numberOfFitParameters();
6057 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6062 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6065 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6075 const PrepRawData *prd{};
6077 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6078 prd = rot->prepRawData();
6088 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6089 if (!splitProb.isSplit()) {
6090 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6094 std::unique_ptr < const RIO_OnTrack > newrot;
6095 double *olderror = state->measurementErrors();
6098 double newerror[5] = {-1,-1,-1,-1,-1};
6099 double newres[2] = {-1,-1};
6101 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6108 newerror[0] = std::sqrt(
covmat(0, 0));
6109 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6110 newerror[1] = std::sqrt(
covmat(1, 1));
6111 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6113 if (
a.cols() != nfitpars) {
6118 for(
int k =0;
k<2;
k++ ){
6119 double oldres =
res[measno+
k];
6120 res[measno+
k] = newres[
k];
6121 err[measno+
k] = newerror[
k];
6123 for (
int i = 0;
i < nfitpars;
i++) {
6124 if (weightderiv(measno+
k,
i) == 0) {
6128 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6130 for (
int j =
i; j < nfitpars; j++) {
6134 weightderiv(measno+
k,
i) *
6135 weightderiv(measno+
k, j) *
6136 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6140 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6144 state->setMeasurement(std::move(newrot));
6145 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.
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 of closest approach of two lines.
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
#define ATH_MSG_WARNING(x)
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
def delta_phi(phi1, phi2)
Gaudi::Property< int > m_maxitPixelROT
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
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.