Loading [MathJax]/extensions/tex2jax.js
 |
ATLAS Offline Software
|
#include <GlobalChi2Fitter.h>
|
| GlobalChi2Fitter (const std::string &, const std::string &, const IInterface *) |
|
virtual | ~ GlobalChi2Fitter () |
|
virtual StatusCode | initialize () override |
|
virtual StatusCode | finalize () override |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const MeasurementSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const PrepRawDataSet &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual std::unique_ptr< Track > | fit (const EventContext &ctx, const Track &, const MeasurementSet &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final |
|
virtual Track * | alignmentFit (AlignmentCache &, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=Trk::nonInteracting) const override |
|
|
Track * | fitIm (const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const |
|
Track * | myfit (const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const |
|
Track * | myfit_helper (Cache &, GXFTrajectory &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const |
|
Track * | mainCombinationStrategy (const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const |
|
Track * | backupCombinationStrategy (const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const |
|
void | makeProtoState (Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const |
|
void | makeProtoStateFromMeasurement (Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const |
|
bool | processTrkVolume (Cache &, const Trk::TrackingVolume *tvol) const |
|
void | addMaterialUpdateTrajectory (Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * >> &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const |
| Given layer information, probe those layers for scatterers and add them to a track. More...
|
|
void | addIDMaterialFast (const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const |
| A faster strategy for adding scatter material to tracks, works only for inner detector tracks. More...
|
|
void | addMaterial (const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const |
|
std::vector< std::unique_ptr< TrackParameters > > | holesearchExtrapolation (const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const |
| Helper method which performs an extrapolation with additional logic for hole search. More...
|
|
std::unique_ptr< const TrackParameters > | makePerigee (Cache &, const TrackParameters &, const ParticleHypothesis) const |
|
std::unique_ptr< const TrackParameters > | makeTrackFindPerigeeParameters (const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
|
std::unique_ptr< GXFTrackState > | makeTrackFindPerigee (const EventContext &, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
|
std::unique_ptr< Track > | makeTrack (const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const |
|
void | fillResidualsAndErrors (const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const |
|
void | tryToConverge (const Cache &cache, GXFTrajectory &trajectory, const int it) const |
|
void | updateSystemWithMaxBremPull (GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const |
|
void | fillDerivatives (GXFTrajectory &traj) const |
|
void | fillFirstLastMeasurement (Cache &cache, GXFTrajectory &trajectory) const |
|
void | fillBfromMeasurements (const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b) const |
|
void | fillAfromMeasurements (const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) const |
|
void | fillAfromScatterers (GXFTrajectory &trajectory, Amg::SymMatrixX &a) const |
|
bool | tryToWeightAfromMaterial (Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2) const |
|
void | compensatePhiWeights (Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) const |
|
FitterStatusCode | runIteration (const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const |
|
FitterStatusCode | updateFitParameters (GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const |
| Method to update peregee parameters, scattering angles, and brems. More...
|
|
void | updatePixelROTs (GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const |
| Update the Pixel ROT using the current trajectory/local track parameters. More...
|
|
GXFTrajectory * | runTrackCleanerSilicon (const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const |
|
void | runTrackCleanerTRT (Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const |
|
PropagationResult | calculateTrackParametersPropagateHelper (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| Helper method that encapsulates calls to the propagator tool in the calculateTrackParameters() method. More...
|
|
PropagationResult | calculateTrackParametersPropagate (const EventContext &, const TrackParameters &, const GXFTrackState &, PropDirection, const MagneticFieldProperties &, bool, bool) const |
| Propagate onto a track state, collecting new track parameters, and optionally the Jacobian and possible holes. More...
|
|
std::vector< std::reference_wrapper< GXFTrackState > > | holeSearchStates (GXFTrajectory &trajectory) const |
| Extracts a collection of track states which are important for hole search. More...
|
|
std::optional< GlobalChi2Fitter::TrackHoleCount > | holeSearchProcess (const EventContext &ctx, const std::vector< std::reference_wrapper< GXFTrackState >> &states) const |
| Conduct a hole search between a list of states, possibly reusing existing information. More...
|
|
void | holeSearchHelper (const std::vector< std::unique_ptr< TrackParameters >> &hc, std::set< Identifier > &id_set, std::set< Identifier > &sct_set, TrackHoleCount &rv, bool count_holes, bool count_dead) const |
| Helper method for the hole search that does the actual counting of holes and dead modules. More...
|
|
FitterStatusCode | calculateTrackParameters (const EventContext &ctx, GXFTrajectory &, bool) const |
|
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > | updateEnergyLoss (const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const |
|
void | calculateTrackErrors (GXFTrajectory &, Amg::SymMatrixX &, bool) const |
|
std::optional< TransportJacobian > | numericalDerivatives (const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const |
|
virtual int | iterationsOfLastFit () const |
|
virtual void | setMinIterations (int) |
|
bool | isMuonTrack (const Track &) const |
|
void | initFieldCache (const EventContext &ctx, Cache &cache) const |
| Initialize a field cache inside a fit cache object. More...
|
|
void | throwFailedToGetTrackingGeomtry () const |
|
bool | ensureValidEntranceCalo (const EventContext &ctx, Cache &cache) const |
|
bool | ensureValidEntranceMuonSpectrometer (const EventContext &ctx, Cache &cache) const |
|
const TrackingGeometry * | trackingGeometry (Cache &cache, const EventContext &ctx) const |
|
const TrackingGeometry * | retrieveTrackingGeometry (const EventContext &ctx) const |
|
|
static std::optional< std::pair< Amg::Vector3D, double > > | addMaterialFindIntersectionDisc (Cache &cache, const DiscSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat) |
| Find the intersection of a set of track parameters onto a disc surface. More...
|
|
static std::optional< std::pair< Amg::Vector3D, double > > | addMaterialFindIntersectionCyl (Cache &cache, const CylinderSurface &surface, const TrackParameters ¶m1, const TrackParameters ¶m2, const ParticleHypothesis mat) |
| Find the intersection of a set of track parameters onto a cylindrical surface. More...
|
|
static void | addMaterialGetLayers (Cache &cache, std::vector< std::pair< const Layer *, const Layer * >> &layers, std::vector< std::pair< const Layer *, const Layer * >> &uplayers, const std::vector< std::unique_ptr< GXFTrackState >> &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat) |
| Collect all possible layers that a given track could have passed through. More...
|
|
static void | makeTrackFillDerivativeMatrix (Cache &, GXFTrajectory &) |
|
static void | calculateDerivatives (GXFTrajectory &) |
|
static bool | correctAngles (double &, double &) |
|
|
ToolHandle< IRIO_OnTrackCreator > | m_ROTcreator {this, "RotCreatorTool", "", ""} |
|
ToolHandle< IRIO_OnTrackCreator > | m_broadROTcreator {this, "BroadRotCreatorTool", "", ""} |
|
ToolHandle< IUpdator > | m_updator {this, "MeasurementUpdateTool", "", ""} |
|
ToolHandle< IExtrapolator > | m_extrapolator {this, "ExtrapolationTool", "Trk::Extrapolator/CosmicsExtrapolator", ""} |
|
ToolHandle< IMultipleScatteringUpdator > | m_scattool {this, "MultipleScatteringTool", "Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator", ""} |
|
ToolHandle< IEnergyLossUpdator > | m_elosstool {this, "EnergyLossTool", "Trk::EnergyLossUpdator/AtlasEnergyLossUpdator", ""} |
|
ToolHandle< IMaterialEffectsUpdator > | m_matupdator {this, "MaterialUpdateTool", "", ""} |
|
ToolHandle< IPropagator > | m_propagator {this, "PropagatorTool", "", ""} |
|
ToolHandle< INavigator > | m_navigator {this, "NavigatorTool", "Trk::Navigator/CosmicsNavigator", ""} |
|
ToolHandle< IResidualPullCalculator > | m_residualPullCalculator {this, "ResidualPullCalculatorTool", "Trk::ResidualPullCalculator/ResidualPullCalculator", ""} |
|
ToolHandle< Trk::ITrkMaterialProviderTool > | m_caloMaterialProvider {this, "CaloMaterialProvider", "Trk::TrkMaterialProviderTool/TrkMaterialProviderTool", ""} |
|
ToolHandle< IMaterialEffectsOnTrackProvider > | m_calotool {this, "MuidTool", "Rec::MuidMaterialEffectsOnTrackProvider/MuidMaterialEffectsOnTrackProvider", ""} |
|
ToolHandle< IMaterialEffectsOnTrackProvider > | m_calotoolparam {this, "MuidToolParam", "", ""} |
|
ToolHandle< IBoundaryCheckTool > | m_boundaryCheckTool {this, "BoundaryCheckTool", "", "Boundary checking tool for detector sensitivities" } |
|
SG::ReadCondHandleKey< TrackingGeometry > | m_trackingGeometryReadKey |
|
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > | m_field_cache_key |
|
const AtlasDetectorID * | m_DetID = nullptr |
|
Gaudi::Property< bool > | m_signedradius {this, "SignedDriftRadius", true} |
|
Gaudi::Property< bool > | m_calomat {this, "MuidMat", false} |
|
Gaudi::Property< bool > | m_extmat {this, "ExtrapolatorMaterial", true} |
|
Gaudi::Property< bool > | m_fillderivmatrix {this, "FillDerivativeMatrix", false} |
|
Gaudi::Property< bool > | m_straightlineprop {this, "StraightLine", true} |
|
Gaudi::Property< bool > | m_extensioncuts {this, "TRTExtensionCuts", true} |
|
Gaudi::Property< bool > | m_sirecal {this, "RecalibrateSilicon", false} |
|
Gaudi::Property< bool > | m_trtrecal {this, "RecalibrateTRT", false} |
|
Gaudi::Property< bool > | m_kinkfinding {this, "KinkFinding", false} |
|
Gaudi::Property< bool > | m_decomposesegments {this, "DecomposeSegments", true} |
|
Gaudi::Property< bool > | m_getmaterialfromtrack {this, "GetMaterialFromTrack", true} |
|
Gaudi::Property< bool > | m_domeastrackpar {this, "MeasuredTrackParameters", true} |
|
Gaudi::Property< bool > | m_storemat {this, "StoreMaterialOnTrack", true} |
|
Gaudi::Property< bool > | m_redoderivs {this, "RecalculateDerivatives", false} |
|
Gaudi::Property< bool > | m_reintoutl {this, "ReintegrateOutliers", false} |
|
Gaudi::Property< bool > | m_acceleration {this, "Acceleration", false} |
|
Gaudi::Property< bool > | m_numderiv {this, "NumericalDerivs", false} |
|
Gaudi::Property< bool > | m_fiteloss {this, "FitEnergyLoss", false} |
|
Gaudi::Property< bool > | m_asymeloss {this, "AsymmetricEnergyLoss", true} |
|
Gaudi::Property< bool > | m_useCaloTG {this, "UseCaloTG", false} |
|
Gaudi::Property< bool > | m_rejectLargeNScat {this, "RejectLargeNScat", false} |
|
Gaudi::Property< bool > | m_createSummary {this, "CreateTrackSummary", true} |
|
Gaudi::Property< bool > | m_holeSearch {this, "DoHoleSearch", false} |
|
Gaudi::Property< double > | m_outlcut {this, "OutlierCut", 5.0} |
|
Gaudi::Property< double > | m_p {this, "Momentum", 0.0} |
|
Gaudi::Property< double > | m_chi2cut {this, "TrackChi2PerNDFCut", 1.e15} |
|
Gaudi::Property< double > | m_scalefactor {this, "TRTTubeHitCut", 2.5} |
|
Gaudi::Property< double > | m_minphfcut {this, "MinPHFCut", 0.} |
|
Gaudi::Property< int > | m_maxoutliers {this, "MaxOutliers", 10} |
|
Gaudi::Property< int > | m_maxit {this, "MaxIterations", 30} |
|
Gaudi::Property< int > | m_miniter {this, "MinimumIterations", 1} |
|
Gaudi::Property< int > | m_fixbrem {this, "FixBrem", -1} |
|
Gaudi::Property< int > | m_maxitPixelROT {this, "IterationsToRebuildPixelRots", 0} |
|
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > | m_clusterSplitProbContainer {this, "ClusterSplitProbabilityName", "",""} |
|
Trk::Volume | m_idVolume |
|
std::array< std::atomic< unsigned int >, S_MAX_VALUE > m_fit_status | ATLAS_THREAD_SAFE = {} |
|
Definition at line 156 of file GlobalChi2Fitter.h.
◆ FitterStatusType
Enumerator |
---|
S_FITS | |
S_SUCCESSFUL_FITS | |
S_MAT_INV_FAIL | |
S_NOT_ENOUGH_MEAS | |
S_PROPAGATION_FAIL | |
S_INVALID_ANGLES | |
S_NOT_CONVERGENT | |
S_HIGH_CHI2 | |
S_LOW_MOMENTUM | |
S_MAX_VALUE | |
Definition at line 176 of file GlobalChi2Fitter.h.
◆ GlobalChi2Fitter()
Trk::GlobalChi2Fitter::GlobalChi2Fitter |
( |
const std::string & |
t, |
|
|
const std::string & |
n, |
|
|
const IInterface * |
p |
|
) |
| |
◆ ~ GlobalChi2Fitter()
◆ addIDMaterialFast()
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
For every track, we need to add its scatterers. That is to say, we need to determine which bits of non-active material the particle in question may have passed through and add them to the track. This is generally an expensive operation, but we can cut some corners if the track only consists of inner detector hits. Specifically, we can exploit the layer structure of the detector to find possible material hits more quickly and efficiently than using the standard material adding algorithm, which is addMaterial.
- Parameters
-
[in,out] | cache | General cache object, as used everywhere. |
[in,out] | trajectory | The current state of the track, respresented in the fitter's internal track representation. States may be added to this. |
[in] | parameters | Starting parameters for the material addition step. |
[in] | part | Standard representation of particle type, used to determine the behaviour of the particle as it traverses materials. |
Definition at line 3465 of file GlobalChi2Fitter.cxx.
3477 if (!caloEntranceIsValid) {
3486 cache.m_negdiscs.empty() &&
3487 cache.m_posdiscs.empty() &&
3488 cache.m_barrelcylinders.empty()
3503 cache.m_fastmat =
false;
3504 addMaterial(ctx, cache, trajectory, refpar2, matEffects);
3520 bool hasmat =
false;
3521 int indexoffset = 0, lastmatindex = 0;
3522 std::vector<std::unique_ptr<GXFTrackState>> & oldstates = trajectory.trackStates();
3524 GXFTrackState *firstsistate =
nullptr;
3525 GXFTrackState *lastsistate =
nullptr;
3536 for (
int i = 0;
i < (
int) oldstates.size();
i++) {
3537 if (oldstates[
i]->materialEffects() !=
nullptr) {
3546 if (firstsistate ==
nullptr) {
3547 if (oldstates[
i]->trackParameters() ==
nullptr) {
3548 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3551 oldstates[
i]->associatedSurface(),
3554 trajectory.m_fieldprop,
3558 if (tmppar ==
nullptr)
return;
3560 oldstates[
i]->setTrackParameters(std::move(tmppar));
3562 firstsistate = oldstates[
i].get();
3564 lastsistate = oldstates[
i].get();
3572 if (lastsistate ==
nullptr) {
3573 throw std::logic_error(
"No track state");
3582 if (lastsistate->trackParameters() ==
nullptr) {
3583 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
3586 lastsistate->associatedSurface(),
3588 trajectory.m_fieldprop,
3592 if (tmppar ==
nullptr)
return;
3594 lastsistate->setTrackParameters(std::move(tmppar));
3603 refpar = lastsistate->trackParameters();
3604 indexoffset = lastmatindex;
3606 refpar = firstsistate->trackParameters();
3619 std::vector<std::pair<const Layer *, const Layer *>>
layers;
3620 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = trajectory.upstreamMaterialLayers();
◆ addMaterial()
Definition at line 3633 of file GlobalChi2Fitter.cxx.
3640 if (refpar2 ==
nullptr) {
3643 const MeasurementBase *firstmuonhit =
nullptr;
3644 const MeasurementBase *lastmuonhit =
nullptr;
3645 const MeasurementBase *firstidhit =
3647 const MeasurementBase *lastidhit =
nullptr;
3648 const MeasurementBase *firsthit =
nullptr;
3649 const MeasurementBase *lasthit =
nullptr;
3650 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3651 std::vector<std::unique_ptr<GXFTrackState>> matstates;
3652 std::unique_ptr< const std::vector < const TrackStateOnSurface *>,
3653 void (*)(
const std::vector<const TrackStateOnSurface *> *) >
3654 matvec(
nullptr,&Trk::GlobalChi2Fitter::Cache::objVectorDeleter<TrackStateOnSurface>);
3655 bool matvec_used=
false;
3656 std::unique_ptr<TrackParameters> startmatpar1;
3657 std::unique_ptr<TrackParameters> startmatpar2;
3668 int npseudomuon1 = 0;
3669 int npseudomuon2 = 0;
3671 for (
auto & state :
states) {
3674 GXFMaterialEffects *meff = state->materialEffects();
3677 if (firstidhit ==
nullptr) {
3686 if (firsthit ==
nullptr) {
3687 firsthit = state->measurement();
3688 if (cache.m_acceleration) {
3689 if (
tp ==
nullptr) {
3693 state->associatedSurface(),
3699 if (
tp ==
nullptr) {
3703 state->setTrackParameters(std::unique_ptr<const TrackParameters>(
tp));
3709 lasthit = state->measurement();
3715 if (firstidhit ==
nullptr) {
3716 firstidhit = state->measurement();
3719 if ((firstidpar ==
nullptr) && (
tp !=
nullptr)) {
3723 lastidhit = state->measurement();
3724 if (
tp !=
nullptr) {
3729 if (firstsiliconpar ==
nullptr) {
3730 firstsiliconpar =
tp;
3732 lastsiliconpar =
tp;
3744 if (firstmuonhit ==
nullptr) {
3745 firstmuonhit = state->measurement();
3746 if (
tp !=
nullptr) {
3750 lastmuonhit = state->measurement();
3751 if (
tp !=
nullptr) {
3757 if (meff->deltaE() == 0) {
3758 if (firstcalopar ==
nullptr) {
3759 firstcalopar = state->trackParameters();
3761 lastcalopar = state->trackParameters();
3763 if (firstmatpar ==
nullptr) {
3764 firstmatpar = state->trackParameters();
3769 std::unique_ptr<TrackParameters> refpar;
3772 if (trajectory.m_straightline &&
m_p != 0) {
3776 refpar = refpar2->associatedSurface().createUniqueTrackParameters(
3777 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3780 if (firstmatpar !=
nullptr) {
3785 if ((startmatpar1 ==
nullptr) || ((firstidhit !=
nullptr) && (firstmuonhit !=
nullptr))) {
3789 double mass = trajectory.mass();
3791 const AmgVector(5) & newpars = startmatpar2->parameters();
3795 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3796 newpars[0], newpars[1], newpars[2], newpars[3],
3797 sign / std::sqrt(oldp * oldp + 2 * 100 *
MeV * std::sqrt(oldp * oldp +
mass *
mass) + 100 *
MeV * 100 *
MeV),
3801 }
else if (trajectory.m_straightline &&
m_p != 0) {
3805 startmatpar1 = startmatpar1->associatedSurface().createUniqueTrackParameters(
3806 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3812 startmatpar2 = startmatpar2->associatedSurface().createUniqueTrackParameters(
3813 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
3817 if ((firstidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3819 DistanceSolution distsol = firstidhit->associatedSurface().straightLineDistanceEstimate(
3821 refpar->momentum().unit()
3826 if (
distance < 0 && distsol.numberOfSolutions() > 0 && !cache.m_acceleration) {
3827 ATH_MSG_DEBUG(
"Obtaining upstream layers from Extrapolator");
3829 const Surface *destsurf = &firstidhit->associatedSurface();
3830 std::unique_ptr<const TrackParameters> tmppar;
3832 if (firstmuonhit !=
nullptr) {
3834 if (caloEntranceIsValid) {
3837 *cache.m_caloEntrance,
3841 if (tmppar !=
nullptr) {
3842 destsurf = &tmppar->associatedSurface();
3847 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
3852 false, matEffects) );
3855 if (matvec && !matvec->empty()) {
3856 for (
int i = (
int)matvec->size() - 1;
i > -1;
i--) {
3857 const MaterialEffectsBase *meb = (*matvec)[
i]->materialEffectsOnTrack();
3860 const MaterialEffectsOnTrack *meot =
static_cast < const MaterialEffectsOnTrack *
>(meb);
3861 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3862 const TrackParameters * newpars = (*matvec)[
i]->trackParameters() !=
nullptr ? (*matvec)[
i]->trackParameters()->clone() :
nullptr;
3863 meff->setSigmaDeltaE(0);
3864 matstates.push_back(std::make_unique<GXFTrackState>(
3866 std::unique_ptr<const TrackParameters>(newpars)
3876 if ((lastidhit !=
nullptr) && trajectory.numberOfSiliconHits() > 0 && cache.m_idmat) {
3877 DistanceSolution distsol = lastidhit->associatedSurface().straightLineDistanceEstimate(
3879 refpar->momentum().unit()
3884 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
3885 ATH_MSG_DEBUG(
"Obtaining downstream ID layers from Extrapolator");
3886 const Surface *destsurf = &lastidhit->associatedSurface();
3887 std::unique_ptr<const TrackParameters> tmppar;
3888 std::unique_ptr<Surface> calosurf;
3889 if (firstmuonhit !=
nullptr) {
3891 if (caloEntranceIsValid) {
3894 *cache.m_caloEntrance,
3899 if (tmppar !=
nullptr) {
3900 const CylinderSurface *cylcalosurf =
nullptr;
3903 cylcalosurf =
static_cast<const CylinderSurface *
>(&tmppar->associatedSurface());
3905 const DiscSurface *disccalosurf =
nullptr;
3908 disccalosurf =
static_cast<const DiscSurface *
>(&tmppar->associatedSurface());
3910 if (cylcalosurf !=
nullptr) {
3912 const CylinderBounds & cylbounds = cylcalosurf->bounds();
3913 double radius = cylbounds.r();
3914 double hlength = cylbounds.halflengthZ();
3915 calosurf = std::make_unique<CylinderSurface>(trans,
radius - 1, hlength);
3916 }
else if (disccalosurf !=
nullptr) {
3918 disccalosurf->center().z() > 0 ?
3919 disccalosurf->center().z() - 1 :
3920 disccalosurf->center().z() + 1
3924 disccalosurf->center().x(),
3925 disccalosurf->center().y(),
3930 trans.translation() << newpos;
3932 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&disccalosurf->bounds());
3933 double rmin = discbounds->rMin();
3934 double rmax = discbounds->rMax();
3935 calosurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
3937 destsurf = calosurf.release();
3941 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
3943 ctx, *startmatpar2, *destsurf,
alongMomentum,
false, matEffects));
3944 matvec_used =
false;
3946 if (matvec && !matvec->empty()) {
3947 for (
const auto &
i : *matvec) {
3952 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
3953 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
3954 if (cache.m_fiteloss && (meot->energyLoss() !=
nullptr)) {
3955 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
3958 if (matEffects ==
electron && cache.m_asymeloss) {
3959 meff->setDeltaE(-5);
3961 if (trajectory.numberOfTRTHits() == 0) {
3962 meff->setScatteringSigmas(0, 0);
3965 meff->setSigmaDeltaE(50);
3968 const TrackParameters * newparams =
i->trackParameters() !=
nullptr ?
i->trackParameters()->clone() :
nullptr;
3970 matstates.push_back(std::make_unique<GXFTrackState>(
3972 std::unique_ptr<const TrackParameters>(newparams)
3984 if (cache.m_calomat && (firstmuonhit !=
nullptr) && (firstidhit !=
nullptr)) {
3987 std::vector<MaterialEffectsOnTrack> calomeots =
m_calotool->extrapolationSurfacesAndEffects(
3991 firstmuonhit->associatedSurface(),
3996 if (calomeots.empty()) {
4001 std::unique_ptr<const TrackParameters> prevtrackpars =
unique_clone(lastidpar);
4002 if (lasthit == lastmuonhit) {
4003 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4006 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4009 calomeots[
i].associatedSurface(),
4012 trajectory.m_fieldprop,
4016 if (layerpar ==
nullptr) {
4021 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4024 lastcalopar = layerpar.get();
4028 double qoverp = layerpar->parameters()[
Trk::qOverP];
4029 double qoverpbrem = 0;
4033 (firstmuonpar !=
nullptr) &&
4034 std::abs(firstmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4036 qoverpbrem = firstmuonpar->parameters()[
Trk::qOverP];
4038 double sign = (qoverp > 0) ? 1 : -1;
4039 qoverpbrem =
sign / (1 / std::abs(qoverp) - std::abs(calomeots[
i].energyLoss()->deltaE()));
4042 const AmgVector(5) & newpar = layerpar->parameters();
4044 layerpar = layerpar->associatedSurface().createUniqueTrackParameters(
4045 newpar[0], newpar[1], newpar[2], newpar[3], qoverpbrem, std::nullopt
4047 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4050 matstates.push_back(std::make_unique<GXFTrackState>(
4052 std::unique_ptr<const TrackParameters>(layerpar !=
nullptr ? layerpar->clone() :
nullptr)
4054 prevtrackpars = std::move(layerpar);
4059 firsthit == firstmuonhit &&
4060 (!cache.m_getmaterialfromtrack || lasthit == lastidhit)
4063 for (
int i = 0;
i < (
int) calomeots.size();
i++) {
4065 std::unique_ptr<const TrackParameters> layerpar(
m_propagator->propagateParameters(
4068 calomeots[
i].associatedSurface(),
4071 trajectory.m_fieldprop,
4075 if (layerpar ==
nullptr) {
4080 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(calomeots[
i]);
4089 double qoverpbrem = layerpar->parameters()[
Trk::qOverP];
4094 (lastmuonpar !=
nullptr) &&
4095 std::abs(lastmuonpar->parameters()[
Trk::qOverP]) > 1.e-9
4099 double sign = (qoverpbrem > 0) ? 1 : -1;
4100 qoverp =
sign / (1 / std::abs(qoverpbrem) + std::abs(calomeots[
i].energyLoss()->deltaE()));
4103 meff->setdelta_p(1000 * (qoverpbrem - qoverp));
4104 const AmgVector(5) & newpar = layerpar->parameters();
4106 prevtrackpars = layerpar->associatedSurface().createUniqueTrackParameters(
4107 newpar[0], newpar[1], newpar[2], newpar[3], qoverp, std::nullopt
4111 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(std::move(meff), std::move(layerpar)));
4116 if (lasthit == lastmuonhit && cache.m_extmat) {
4117 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4119 if (lastcalopar !=
nullptr) {
4121 if (msEntranceIsValid) {
4122 if (cache.m_msEntrance->inside(lastcalopar->position())) {
4125 *cache.m_msEntrance,
4129 if (muonpar1 !=
nullptr) {
4137 rot.col(2) = trackdir;
4139 trans.linear().matrix() << rot;
4140 trans.translation() << muonpar1->
position() - .1 * trackdir;
4141 PlaneSurface curvlinsurf(trans);
4143 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4151 if (curvlinpar !=
nullptr) {
4152 muonpar1 = std::move(curvlinpar);
4156 muonpar1 = std::unique_ptr<TrackParameters>(lastcalopar->clone());
4160 muonpar1 = std::unique_ptr<TrackParameters>(refpar->clone());
4163 DistanceSolution distsol;
4165 if (muonpar1 !=
nullptr) {
4166 distsol = lastmuonhit->associatedSurface().straightLineDistanceEstimate(
4174 if ((
distance > 0) and(distsol.numberOfSolutions() >
4175 0) and (firstmuonhit !=
nullptr)) {
4176 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4183 if (distsol.numberOfSolutions() == 1) {
4185 }
else if (distsol.numberOfSolutions() == 2) {
4187 std::abs(distsol.first()) < std::abs(distsol.second()) ?
4193 if (
distance < 0 && distsol.numberOfSolutions() > 0 && (firstidhit ==
nullptr)) {
4194 if (firstmuonpar !=
nullptr) {
4197 if (trajectory.m_straightline &&
m_p != 0) {
4202 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4205 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4208 firstmuonhit->associatedSurface(),
4211 trajectory.m_fieldprop,
4215 if (tmppar !=
nullptr) {
4216 muonpar1 = std::move(tmppar);
4222 ATH_MSG_DEBUG(
"Obtaining downstream layers from Extrapolator");
4223 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4226 states.back()->associatedSurface(),
4230 matvec_used =
false;
4237 if (matvec && !matvec->empty()) {
4238 for (
int j = 0; j < (
int) matvec->size(); j++) {
4239 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4243 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4244 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4247 !trajectory.m_straightline &&
4248 (meot->energyLoss() !=
nullptr) &&
4249 std::abs(meff->deltaE()) > 25 &&
4250 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4252 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4255 const TrackParameters * newparams = (*matvec)[j]->trackParameters() !=
nullptr ? (*matvec)[j]->trackParameters()->clone() :
nullptr;
4257 matstates.push_back(std::make_unique<GXFTrackState>(
4259 std::unique_ptr<const TrackParameters>(newparams)
4269 if (firsthit == firstmuonhit && cache.m_extmat && (firstcalopar !=
nullptr)) {
4270 std::unique_ptr<const Trk::TrackParameters> muonpar1;
4273 if (msEntranceIsValid) {
4274 if (cache.m_msEntrance->inside(firstcalopar->position())) {
4277 *cache.m_msEntrance,
4281 if (muonpar1 !=
nullptr) {
4289 rot.col(2) = trackdir;
4291 trans.linear().matrix() << rot;
4292 trans.translation() << muonpar1->
position() - .1 * trackdir;
4293 PlaneSurface curvlinsurf(trans);
4295 std::unique_ptr<const TrackParameters> curvlinpar(
m_extrapolator->extrapolateDirectly(
4303 if (curvlinpar !=
nullptr) {
4304 muonpar1 = std::move(curvlinpar);
4308 muonpar1 = std::unique_ptr<const TrackParameters>(firstcalopar->clone());
4312 DistanceSolution distsol;
4314 if (muonpar1 !=
nullptr) {
4315 distsol = firstmuonhit->associatedSurface().straightLineDistanceEstimate(
4323 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4325 ATH_MSG_DEBUG(
"Collecting upstream muon material from extrapolator");
4326 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
4329 states[0]->associatedSurface(),
4333 matvec_used =
false;
4335 if (matvec && !matvec->empty()) {
4336 ATH_MSG_DEBUG(
"Retrieved " << matvec->size() <<
" material states");
4338 for (
int j = 0; j < (
int) matvec->size(); j++) {
4339 const MaterialEffectsBase *meb = (*matvec)[j]->materialEffectsOnTrack();
4341 if (meb !=
nullptr) {
4345 const MaterialEffectsOnTrack *meot =
static_cast<const MaterialEffectsOnTrack *
>(meb);
4346 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>(*meot);
4349 !trajectory.m_straightline &&
4350 (meot->energyLoss() !=
nullptr) &&
4351 std::abs(meff->deltaE()) > 25 &&
4352 std::abs((*matvec)[j]->trackParameters()->position().z()) < 13000
4354 meff->setSigmaDeltaE(meot->energyLoss()->sigmaDeltaE());
4358 (*matvec)[j]->trackParameters() !=
nullptr
4359 ? (*matvec)[j]->trackParameters()->clone()
4362 matstates.insert(matstates.begin(), std::make_unique<GXFTrackState>(
4364 std::unique_ptr<const TrackParameters>(tmpparams)
4377 std::vector<std::unique_ptr<GXFTrackState>> & newstates =
states;
4378 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(newstates);
4381 newstates.reserve(oldstates.size() + matstates.size());
4384 int firstlayerno = -1;
4386 if (cache.m_acceleration) {
4387 trajectory.addBasicState(std::move(oldstates[0]));
4393 for (
int i = cache.m_acceleration ? 1 : 0;
i < (
int) oldstates.size();
i++) {
4394 bool addlayer =
true;
4396 while (addlayer && layerno < (
int) matstates.size()) {
4398 const TrackParameters *layerpar = matstates[layerno]->trackParameters();
4400 DistanceSolution distsol = oldstates[
i]->associatedSurface().straightLineDistanceEstimate(
4401 layerpar->position(),
4402 layerpar->momentum().unit()
4407 if (
distance > 0 && distsol.numberOfSolutions() > 0) {
4412 double cylinderradius = layerpar->associatedSurface().bounds().r();
4413 double trackimpact = std::abs(-refpar->position().x() * sinphi + refpar->position().y() * cosphi);
4415 if (trackimpact > cylinderradius - 5 *
mm) {
4421 if (
i == (
int) oldstates.size() - 1) {
4426 GXFMaterialEffects *meff = matstates[layerno]->materialEffects();
4428 if (meff->sigmaDeltaPhi() > .4 || (meff->sigmaDeltaPhi() == 0 && meff->sigmaDeltaE() <= 0)) {
4429 if (meff->sigmaDeltaPhi() > .4) {
4430 ATH_MSG_DEBUG(
"Material state with excessive scattering, skipping it");
4433 if (meff->sigmaDeltaPhi() == 0) {
4441 if (firstlayerno < 0) {
4442 firstlayerno = layerno;
4445 trajectory.addMaterialState(std::move(matstates[layerno]));
4447 if ((layerpar !=
nullptr) && matEffects !=
pion && matEffects !=
muon) {
4448 const TrackingVolume *tvol =
m_navigator->volume(ctx,layerpar->position());
4449 const Layer *lay =
nullptr;
4451 if (tvol !=
nullptr) {
4452 lay = (tvol->closestMaterialLayer(layerpar->position(),layerpar->momentum().normalized())).
object;
4455 const MaterialProperties *matprop = lay !=
nullptr ? lay->fullUpdateMaterialProperties(*layerpar) :
nullptr;
4456 meff->setMaterialProperties(matprop);
4462 trajectory.addBasicState(std::move(oldstates[
i]));
4465 ATH_MSG_DEBUG(
"Total X0: " << trajectory.totalX0() <<
" total eloss: " << trajectory.totalEnergyLoss());
4467 if (matvec_used) cache.m_matTempStore.push_back( std::move(matvec) );
◆ addMaterialFindIntersectionCyl()
Find the intersection of a set of track parameters onto a cylindrical surface.
See addMaterialFindIntersectionDisc for more information.
- Note
- This method can probably be replaced entirely by the straight line intersection method of the appropriate Surface subclass.
Definition at line 2995 of file GlobalChi2Fitter.cxx.
3009 const double *
pos = parforextrap.position().data();
3011 cache.m_field_cache.getFieldZR(
pos,
field);
3017 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
3018 double xc = parforextrap.position().x() -
r * sinphi;
3019 double yc = parforextrap.position().y() +
r * cosphi;
3020 double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
3021 double z0 = parforextrap.position().z();
3022 double d = xc * xc + yc * yc;
3023 double rcyl = surf.bounds().r();
3024 double mysqrt = ((
r + rcyl) * (
r + rcyl) -
d) * (
d - (
r - rcyl) * (
r - rcyl));
3030 mysqrt = std::sqrt(mysqrt);
3031 double firstterm = xc / 2 + (xc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3032 double secondterm = (mysqrt * yc) / (2 *
d);
3033 double x1 = firstterm + secondterm;
3034 double x2 = firstterm - secondterm;
3035 firstterm = yc / 2 + (yc * (rcyl * rcyl -
r *
r)) / (2 *
d);
3036 secondterm = (mysqrt * xc) / (2 *
d);
3037 double y1 = firstterm - secondterm;
3038 double y2 = firstterm + secondterm;
3039 double x = parforextrap.position().x();
3040 double y = parforextrap.position().y();
3041 double dist1 = (
x -
x1) * (
x -
x1) + (
y -
y1) * (
y -
y1);
3042 double dist2 = (
x -
x2) * (
x -
x2) + (
y -
y2) * (
y -
y2);
3044 if (dist1 < dist2) {
3052 double phi1 = std::atan2(
y - yc,
x - xc);
3055 double delta_z =
r * deltaphi / tantheta;
3056 double z =
z0 + delta_z;
3060 if (std::abs(
z - surf.center().z()) > surf.bounds().halflengthZ()) {
3069 double costracksurf = std::abs(normal.unit().dot(trackdir));
3071 return std::make_pair(
intersect, costracksurf);
◆ addMaterialFindIntersectionDisc()
Find the intersection of a set of track parameters onto a disc surface.
Calculates the intersection from a point and momentum in space onto a disc surface which represents a disc-shaped layer in the detector. The position of the intersection can be used to find materials in that layer at that position.
- Parameters
-
[in] | cache | The standard GX2F cache. |
[in] | surface | The surface to intersect with. |
[in] | param1 | The main track parameters to calculate the intersection from. |
[in] | param2 | A secondary set of parameters used for electrons. The purpose of this is not known to us at this time. |
[in] | mat | A particle hypothesis describing the behaviour of the particle. |
- Returns
- Nothing if the intersection failed (i.e. there was no intersection), otherwise both an intersection positition as well as the angle of inflection.
- Note
- This method can probably be replaced entirely by the straight line intersection method of the appropriate Surface subclass.
Definition at line 2950 of file GlobalChi2Fitter.cxx.
2962 const double *
pos = parforextrap.position().data();
2964 cache.m_field_cache.getFieldZR(
pos,
field);
2972 double r = (std::abs(currentqoverp) > 1
e-10) ? -sintheta / (currentqoverp * 300. *
field[2]) : 1e6;
2973 double xc = parforextrap.position().x() -
r * sinphi;
2974 double yc = parforextrap.position().y() +
r * cosphi;
2975 double phi0 = std::atan2(parforextrap.position().y() - yc, parforextrap.position().x() - xc);
2976 double z0 = parforextrap.position().z();
2977 double delta_s = (surf.center().z() -
z0) / costheta;
2983 const DiscBounds *discbounds = (
const DiscBounds *) (&surf.bounds());
2985 if (
perp > discbounds->rMax() ||
perp < discbounds->rMin()) {
2989 double costracksurf = std::abs(costheta);
2991 return std::make_pair(
intersect, costracksurf);
◆ addMaterialGetLayers()
Collect all possible layers that a given track could have passed through.
If we are to use layer information to determine possible scatterer hits, we must first gather those layers. That's what this method does. It looks for disc and barrel cylinder layers that the given track might have crossed and collects them into output vectors. One contains layers between states on the track, and the upstream layers lie before the first state of the track.
- Parameters
-
[in,out] | cache | General cache object. |
[out] | layers | Output vector for layers. |
[out] | uplayers | Output vector for upstream layers, which lie before the first hit in the track. |
[in] | states | A list of track states on the track. |
[in] | first | The first track state. |
[in] | last | The last track state. |
[in] | refpar | Reference parameters from which to extrapolate. |
[in] | hasmat | Are there any existing materials on this track? |
Definition at line 3285 of file GlobalChi2Fitter.cxx.
3298 upstreamlayers.reserve(5);
3305 double firstz = firstsistate.trackParameters()->position().z();
3306 double firstr = firstsistate.trackParameters()->position().perp();
3307 double firstz2 = hasmat ? lastsistate.trackParameters()->position().z() : firstsistate.trackParameters()->position().z();
3308 double firstr2 = hasmat ? lastsistate.trackParameters()->position().perp() : firstsistate.trackParameters()->position().perp();
3310 GXFTrackState *firststate = oldstates.front().get();
3311 GXFTrackState *laststate = oldstates.back().get();
3317 double lastz = laststate->position().z();
3318 double lastr = laststate->position().perp();
3320 const Layer *startlayer = firststate->associatedSurface().associatedLayer();
3321 const Layer *startlayer2 = hasmat ? lastsistate.associatedSurface().associatedLayer() :
nullptr;
3322 const Layer *endlayer = laststate->associatedSurface().associatedLayer();
3325 double slope = (tantheta != 0) ? 1 / tantheta : 0;
3331 std::vector < const Layer *>::const_iterator
it;
3332 std::vector < const Layer *>::const_iterator itend;
3340 it = cache.m_posdiscs.begin();
3341 itend = cache.m_posdiscs.end();
3343 it = cache.m_negdiscs.begin();
3344 itend = cache.m_negdiscs.end();
3350 for (;
it != itend; ++
it) {
3355 if (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(lastz)) {
3363 const DiscBounds *discbounds = (
const DiscBounds *) (&(*it)->surfaceRepresentation().bounds());
3368 if (discbounds->rMax() < firstr || discbounds->rMin() > lastr) {
3372 double rintersect = firstr + ((*it)->surfaceRepresentation().center().z() - firstz) / slope;
3375 rintersect < discbounds->rMin() - 50 ||
3376 rintersect > discbounds->rMax() + 50
3386 if ((*
it) == endlayer) {
3398 std::abs((*it)->surfaceRepresentation().center().z()) < std::abs(firstz) ||
3401 upstreamlayers.emplace_back((
Layer *)
nullptr, (*it));
3408 (*
it) != startlayer &&
3409 (std::abs((*it)->surfaceRepresentation().center().z()) > std::abs(firstz2) ||
3410 (*it) == startlayer2)
3420 for (
const auto *barrelcylinder : cache.m_barrelcylinders) {
3424 if (barrelcylinder->surfaceRepresentation().bounds().r() > lastr) {
3431 double zintersect = firstz + (barrelcylinder->surfaceRepresentation().bounds().r() - firstr) * slope;
3433 if (std::abs(zintersect - barrelcylinder->surfaceRepresentation().center().z()) >
3434 ((
const CylinderSurface*)(&barrelcylinder->surfaceRepresentation()))->bounds().halflengthZ() + 50) {
3438 if (barrelcylinder == endlayer) {
3445 if (barrelcylinder->surfaceRepresentation().bounds().r() < firstr ||
3446 barrelcylinder == startlayer) {
3447 upstreamlayers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
3450 if (barrelcylinder != startlayer &&
3451 (barrelcylinder->surfaceRepresentation().bounds().r() > firstr2 ||
3452 barrelcylinder == startlayer2)) {
3453 layers.emplace_back(barrelcylinder, (
Layer*)
nullptr);
◆ addMaterialUpdateTrajectory()
Given layer information, probe those layers for scatterers and add them to a track.
This is the meat of the pudding, if you will. Given the information that we have about layers, go through them all and find any possible material hits that we need to add to the track.
- Parameters
-
[in,out] | cache | General cache object. |
[in,out] | track | The track object as it exists now in IR. |
[in] | offset | The first state after any existing materials. |
[in] | layers | The list of layers. |
[in] | ref1 | The first set of reference parameters. |
[in] | ref2 | The second set of reference parameters. |
[in] | mat | The particle hypothesis describing the track behaviour. |
- Note
- Attentive readers may wonder why we pass this function a vector of layers, but not a vector of upstream layers. The reason for this is that the vector of upstream layers is also a member of the cache object.
Definition at line 3074 of file GlobalChi2Fitter.cxx.
3083 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
3084 std::vector<std::unique_ptr<GXFTrackState>> oldstates = std::move(
states);
3095 for (
int i = 0;
i <= indexoffset;
i++) {
3096 trajectory.addBasicState(std::move(oldstates[
i]));
3105 for (
int i = indexoffset + 1;
i < (
int) oldstates.size();
i++) {
3106 double rmeas = oldstates[
i]->position().perp();
3107 double zmeas = oldstates[
i]->position().z();
3115 while (layerindex < (
int)
layers.size()) {
3117 double costracksurf = 0.0;
3132 const CylinderSurface *cylsurf = (
const CylinderSurface *) (&
layer->surfaceRepresentation());
3138 if (oldstates[
i]->trackParameters() !=
nullptr) {
3139 double rlayer = cylsurf->bounds().r();
3140 if (std::abs(rmeas - rlayer) < std::abs(parforextrap->position().perp() - rlayer)) {
3141 parforextrap = oldstates[
i]->trackParameters();
3157 if (cylsurf->bounds().r() > rmeas)
break;
3165 const DiscSurface *discsurf = (
const DiscSurface *) (&
layer->surfaceRepresentation());
3167 if (oldstates[
i]->trackParameters() !=
nullptr) {
3168 double zlayer = discsurf->center().z();
3169 if (std::abs(zmeas - zlayer) < std::abs(parforextrap->position().z() - zlayer)) {
3170 parforextrap = oldstates[
i]->trackParameters();
3181 if (std::abs(discsurf->center().z()) > std::abs(zmeas))
break;
3183 throw std::logic_error(
"Unhandled surface.");
3190 const MaterialProperties *matprop =
layer->layerMaterialProperties()->fullMaterial(
intersect);
3191 if (matprop ==
nullptr) {
3200 double X0 = matprop->thicknessInX0();
3204 double actualx0 =
X0 / costracksurf;
3205 double de = -std::abs(
3206 (matprop->thickness() / costracksurf) *
3209 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3212 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
3214 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3218 std::unique_ptr<GXFMaterialEffects> meff = std::make_unique<GXFMaterialEffects>();
3219 meff->setDeltaE(de);
3220 meff->setScatteringSigmas(sigmascat / sintheta, sigmascat);
3221 meff->setX0(actualx0);
3222 meff->setSurface(&
layer->surfaceRepresentation());
3223 meff->setMaterialProperties(matprop);
3229 std::unique_ptr<EnergyLoss> eloss;
3231 if (cache.m_fiteloss || (matEffects ==
electron && cache.m_asymeloss)) {
3232 eloss = std::make_unique<EnergyLoss>(
m_elosstool->energyLoss(
3234 (
m_p != 0.0 ? std::abs(
m_p) : std::abs(1. / currentqoverp)),
3239 if (eloss !=
nullptr) {
3240 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
3244 if (matEffects ==
electron && cache.m_asymeloss) {
3245 meff->setDeltaE(-5);
3246 if (trajectory.numberOfTRTHits() == 0) {
3247 meff->setScatteringSigmas(0, 0);
3250 meff->setSigmaDeltaE(50);
3251 if (eloss !=
nullptr) {
3252 meff->setSigmaDeltaEPos(eloss->sigmaPlusDeltaE());
3253 meff->setSigmaDeltaENeg(eloss->sigmaMinusDeltaE());
3258 "X0: " << meff->x0() <<
" qoverp: " << currentqoverp <<
3259 " sigmascat " << meff->sigmaDeltaTheta() <<
" eloss: " << meff->deltaE() <<
3260 " sigma eloss: " << meff->sigmaDeltaE()
3267 std::unique_ptr<GXFTrackState> matstate = std::make_unique<GXFTrackState>(
3269 std::unique_ptr<const TrackParameters>()
3272 trajectory.addMaterialState(std::move(matstate));
3281 trajectory.addBasicState(std::move(oldstates[
i]));
◆ alignmentFit()
Definition at line 1819 of file GlobalChi2Fitter.cxx.
1825 const EventContext& ctx = Gaudi::Hive::currentContext();
1829 delete alignCache.m_derivMatrix;
1830 alignCache.m_derivMatrix =
nullptr;
1832 delete alignCache.m_fullCovarianceMatrix;
1833 alignCache.m_fullCovarianceMatrix =
nullptr;
1834 alignCache.m_iterationsOfLastFit = 0;
1837 fitIm(ctx, cache, inputTrack, runOutlier, matEffects);
1838 if(newTrack !=
nullptr){
1839 if(cache.m_derivmat.size() != 0)
1840 alignCache.m_derivMatrix =
new Amg::MatrixX(cache.m_derivmat);
1841 if(cache.m_fullcovmat.size() != 0)
1842 alignCache.m_fullCovarianceMatrix =
new Amg::MatrixX(cache.m_fullcovmat);
1843 alignCache.m_iterationsOfLastFit = cache.m_lastiter;
◆ backupCombinationStrategy()
Definition at line 1242 of file GlobalChi2Fitter.cxx.
1250 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::backupCombinationStrategy");
1253 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
1263 const auto *
const pParametersVector = indettrack->trackParameters();
1266 if (pParametersVector->size() > 1)
1267 firstidpar = (*pParametersVector)[1];
1269 firstidpar = pParametersVector->back();
1271 std::unique_ptr<const TrackParameters> lastidpar =
nullptr;
1272 if ((firstidpar !=
nullptr) && (cache.m_caloEntrance !=
nullptr))
1276 if (lastidpar ==
nullptr) {
1277 lastidpar.reset(pParametersVector->back()->clone());
1280 std::unique_ptr < const TrackParameters > firstscatpar;
1281 std::unique_ptr < const TrackParameters > lastscatpar;
1282 std::unique_ptr < const TrackParameters > elosspar;
1287 lastidpar->position(),
1288 lastidpar->momentum(),
1290 lastidpar->position()
1297 calomeots[0].associatedSurface(),
1300 trajectory.m_fieldprop,
1303 if (!firstscatpar) {
1307 std::unique_ptr<const TrackParameters> tmppar(
1311 calomeots[1].associatedSurface(),
1314 trajectory.m_fieldprop,
1324 double oldp = std::abs(1 / tmppar->parameters()[
Trk::qOverP]);
1325 double de = std::abs(calomeots[1].energyLoss()->deltaE());
1327 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
1333 double newqoverp =
sign / std::sqrt(newp2);
1338 tmppar->associatedSurface().createUniqueTrackParameters(
1345 calomeots[2].associatedSurface(),
1348 trajectory.m_fieldprop,
1358 calomeots[2].associatedSurface(),
1361 trajectory.m_fieldprop,
1372 calomeots[1].associatedSurface(),
1375 trajectory.m_fieldprop,
1384 double newqoverp =
sign /
1385 (1. / std::abs(elosspar->parameters()[
Trk::qOverP]) +
1386 std::abs(calomeots[1].energyLoss()->deltaE()));
1390 std::unique_ptr<const TrackParameters>tmppar(
1391 elosspar->associatedSurface().createUniqueTrackParameters(
1399 calomeots[0].associatedSurface(),
1402 trajectory.m_fieldprop,
1406 if (!firstscatpar) {
1411 for (; itStates != endState; ++itStates) {
1419 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
1421 cache.m_idmat =
false;
1432 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1433 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1434 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1438 sigmadp = calomeots[1].energyLoss()->sigmaDeltaE();
1439 elossmeff->setSigmaDeltaE(sigmadp);
1442 elossmeff->setdelta_p(
dp);
1444 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1445 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1446 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1448 GXFTrackState *secondscatstate = trajectory.trackStates().back().get();
1449 const Surface *triggersurf1 =
nullptr;
1450 const Surface *triggersurf2 =
nullptr;
1454 bool seenmdt =
false;
1455 bool mdtbetweenphihits =
false;
1459 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1460 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1461 (!firstismuon ? ++itStates2 : --itStates2)
1464 ((*itStates2)->measurementOnTrack() ==
nullptr) ||
1469 const auto *
const pMeasurement = (*itStates2)->measurementOnTrack();
1470 const Surface *surf = &pMeasurement->associatedSurface();
1474 if (isCompetingRIOsOnTrack) {
1476 rot = &crot->rioOnTrack(0);
1479 rot =
static_cast<const RIO_OnTrack *
>(pMeasurement);
1482 if ((rot !=
nullptr) &&
m_DetID->
is_mdt(rot->identify()) && (triggersurf1 !=
nullptr)) {
1486 (rot !=
nullptr) && (
1492 Amg::Vector3D measdir = surf->transform().rotation().col(0);
1494 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1496 bool measphi = std::abs(dotprod1) <= .5 && std::abs(dotprod2) <= .5;
1500 (*itStates2)->trackParameters() !=
nullptr ?
1501 (*itStates2)->trackParameters()->position() :
1502 rot->globalPosition();
1503 if (triggersurf1 !=
nullptr) {
1504 triggerpos2 = thispos;
1505 triggersurf2 = surf;
1507 mdtbetweenphihits =
true;
1510 triggerpos1 = thispos;
1511 triggersurf1 = surf;
1517 double mdttrig1 = 999999;
1518 double mdttrig2 = 999999;
1519 const Surface *mdtsurf1 =
nullptr;
1520 const Surface *mdtsurf2 =
nullptr;
1523 itStates2 = (!firstismuon ? beginStates2 : endState - 1);
1524 itStates2 != (!firstismuon ? endState2 : beginStates - 1);
1525 (!firstismuon ? ++itStates2 : --itStates2)
1527 const Surface *surf =
nullptr;
1529 ((*itStates2)->measurementOnTrack() !=
nullptr) &&
1532 surf = &(*itStates2)->measurementOnTrack()->associatedSurface();
1535 if (surf ==
nullptr) {
1538 const auto *
const pThisMeasurement = (*itStates2)->measurementOnTrack();
1542 if (isCompetingRioOnTrack) {
1544 rot = &crot->rioOnTrack(0);
1547 rot =
static_cast<const RIO_OnTrack *
>(pThisMeasurement);
1550 const bool thisismdt = rot and
m_DetID->
is_mdt(rot->identify());
1553 (*itStates2)->trackParameters() !=
nullptr ?
1554 (*itStates2)->trackParameters()->position() :
1555 pThisMeasurement->globalPosition();
1556 if (triggerpos1.mag() > 1 && (globpos - triggerpos1).
mag() < mdttrig1) {
1557 mdttrig1 = (globpos - triggerpos1).
mag();
1560 if (triggerpos2.mag() > 1 && (globpos - triggerpos2).
mag() < mdttrig2) {
1561 mdttrig2 = (globpos - triggerpos2).
mag();
1567 GXFTrackState * firstpseudostate =
nullptr;
1568 std::vector<GXFTrackState *> outlierstates;
1569 std::vector<GXFTrackState *> outlierstates2;
1571 outlierstates.reserve(10);
1572 outlierstates2.reserve(10);
1574 std::unique_ptr<PseudoMeasurementOnTrack> newpseudo;
1576 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1577 const auto *
const pMeasurement{(*itStates2)->measurementOnTrack()};
1579 bool isStraightLine =
1580 pMeasurement !=
nullptr ?
1587 (newpseudo ==
nullptr) && (
1588 (itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1589 std::abs(pMeasurement->globalPosition().z()) < 10000
1592 std::unique_ptr<const TrackParameters> par2;
1593 if (((*itStates2)->trackParameters() !=
nullptr) && nphi > 99) {
1594 par2.reset((*itStates2)->trackParameters()->clone());
1598 *secondscatstate->trackParameters(),
1599 pMeasurement->associatedSurface(),
1602 trajectory.m_fieldprop,
1606 if (par2 ==
nullptr) {
1612 newpseudo = std::make_unique<PseudoMeasurementOnTrack>(
1615 par2->associatedSurface()
1618 std::unique_ptr<GXFTrackState> firstpseudo = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(par2));
1625 firstpseudo->setMeasurementErrors(
errors);
1626 firstpseudostate = firstpseudo.get();
1627 trajectory.addMeasurementState(std::move(firstpseudo));
1632 if (isPseudo && !firstismuon) {
1636 if ((**itStates2).materialEffectsOnTrack() !=
nullptr) {
1638 cache.m_idmat =
false;
1646 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1647 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf1 &&
1648 (mdtsurf1 !=
nullptr)
1650 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf1->transform());
1652 transf->translation() << triggerpos1;
1653 StraightLineSurface slsurf(*transf);
1657 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1661 std::unique_ptr<GXFTrackState> pseudostate1 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1668 pseudostate1->setMeasurementErrors(
errors);
1669 outlierstates2.push_back(pseudostate1.get());
1670 trajectory.addMeasurementState(std::move(pseudostate1));
1674 ((**itStates2).measurementOnTrack() !=
nullptr) &&
1675 &(**itStates2).measurementOnTrack()->associatedSurface() == triggersurf2 &&
1676 mdtbetweenphihits &&
1677 (mdtsurf2 !=
nullptr)
1679 std::unique_ptr<Amg::Transform3D> transf = std::make_unique<Amg::Transform3D>(mdtsurf2->transform());
1680 transf->translation() << triggerpos2;
1681 StraightLineSurface slsurf(*transf);
1685 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1689 std::unique_ptr<GXFTrackState> pseudostate2 = std::make_unique<GXFTrackState>(std::move(newpseudo),
nullptr);
1696 pseudostate2->setMeasurementErrors(
errors);
1698 outlierstates2.push_back(pseudostate2.get());
1699 trajectory.addMeasurementState(std::move(pseudostate2));
1706 trajectory.trackStates().back()->measurementType() ==
TrackState::TGC ||
1708 trajectory.trackStates().back()->measurementType() ==
TrackState::RPC &&
1709 trajectory.trackStates().back()->measuresPhi()
1714 outlierstates.push_back(trajectory.trackStates().back().get());
1715 trajectory.setOutlier((
int) trajectory.trackStates().size() - 1,
true);
1720 trajectory.setNumberOfPerigeeParameters(0);
1724 trajectory.setPrefit(2);
1726 cache.m_matfilled =
true;
1727 bool tmpacc = cache.m_acceleration;
1728 cache.m_acceleration =
false;
1730 std::unique_ptr<Trk::Track> tmp_track(
1731 myfit(ctx, cache, trajectory, *startpar2,
false,
muon));
1732 cache.m_acceleration = tmpacc;
1734 cache.m_matfilled =
false;
1737 trajectory.converged() &&
1738 std::abs(trajectory.residuals().tail<1>()(0) / trajectory.errors().tail<1>()(0)) > 10
1743 if (trajectory.converged()) {
1744 if (firstpseudostate !=
nullptr) {
1749 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1752 par2->associatedSurface()
1754 firstpseudostate->setMeasurement(std::move(newpseudo));
1755 firstpseudostate->setRecalibrated(
false);
1758 for (
int j = 0; j < (
int) trajectory.trackStates().size(); j++) {
1759 for (
const auto &
i : outlierstates2) {
1760 if (trajectory.trackStates()[j].get() ==
i) {
1761 trajectory.setOutlier(j,
true);
1765 for (
const auto &
i : outlierstates) {
1766 if (trajectory.trackStates()[j].get() ==
i) {
1767 trajectory.setOutlier(j,
false);
1773 itStates = (firstismuon ? beginStates2 : endState - 1);
1774 itStates != (firstismuon ? endState2 : beginStates - 1);
1775 (firstismuon ? ++itStates : --itStates)
1781 makeProtoState(cache, trajectory, *itStates, (firstismuon ? -1 : 0));
1785 trajectory.setPrefit(0);
1786 trajectory.setNumberOfPerigeeParameters(5);
1788 cache.m_matfilled =
false;
◆ calculateDerivatives()
void Trk::GlobalChi2Fitter::calculateDerivatives |
( |
GXFTrajectory & |
trajectory | ) |
|
|
staticprivate |
Definition at line 8007 of file GlobalChi2Fitter.cxx.
8008 int nstatesupstream = trajectory.numberOfUpstreamStates();
8009 int nscatupstream = trajectory.numberOfUpstreamScatterers();
8010 int nbremupstream = trajectory.numberOfUpstreamBrems();
8011 int nscats = trajectory.numberOfScatterers();
8012 int nperpars = trajectory.numberOfPerigeeParameters();
8013 int nfitpars = trajectory.numberOfFitParameters();
8015 using Matrix55 = Eigen::Matrix<double, 5, 5>;
8017 Matrix55 initialjac;
8018 initialjac.setZero();
8019 initialjac(4, 4) = 1;
8021 Matrix55 jacvertex(initialjac);
8023 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacscat(trajectory.numberOfScatterers(), initialjac);
8024 std::vector<Matrix55, Eigen::aligned_allocator<Matrix55>> jacbrem(trajectory.numberOfBrems(), initialjac);
8026 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
8027 GXFTrackState *prevstate =
nullptr, *state =
nullptr;
8029 int hit_begin = 0, hit_end = 0, scatno = 0, bremno = 0;
8031 for (
bool forward : {
false,
true}) {
8033 hit_begin = nstatesupstream;
8035 scatno = nscatupstream;
8036 bremno = nbremupstream;
8038 hit_begin = nstatesupstream - 1;
8040 scatno = trajectory.numberOfUpstreamScatterers() - 1;
8041 bremno = trajectory.numberOfUpstreamBrems() - 1;
8045 int hitno = hit_begin;
8046 forward ? (hitno < hit_end) : (hitno >= hit_end);
8047 hitno += (forward ? 1 : -1)
8050 state =
states[hitno].get();
8054 if (fillderivmat && state->derivatives().cols() != nfitpars) {
8055 state->derivatives().resize(5, nfitpars);
8056 state->derivatives().setZero();
8059 int jminscat = 0, jmaxscat = 4, jminbrem = 0, jmaxbrem = 4;
8061 if (hitno == (forward ? hit_end - 1 : 0)) {
8062 if (!fillderivmat) {
8070 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8072 if (hitno == nstatesupstream + (forward ? 0 : -1)) {
8073 jacvertex.block<4, 5>(0, 0) = jac.block<4, 5>(0, 0);
8074 jacvertex(4, 4) = jac(4, 4);
8076 int jmin = 0, jmax = 0, jcnt = 0;
8077 int lp_bgn = 0, lp_end = 0;
8081 jcnt = jmax - jmin + 1;
8083 lp_bgn = forward ? nscatupstream : nscatupstream - 1;
8086 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8088 i == scatno + (forward ? -1 : 1) &&
8089 prevstate !=
nullptr &&
8091 (!trajectory.prefit() || prevstate->materialEffects()->deltaE() == 0)
8093 jacscat[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8094 jacscat[
i](4, 4) = jac(4, 4);
8096 calculateJac(jac, jacscat[
i], jmin, jmax);
8100 Eigen::MatrixXd & derivmat = state->derivatives();
8101 int scatterPos = nperpars + 2 *
i;
8103 derivmat.block<4, 2>(0, scatterPos) = (forward ? 1 : -1) * jacscat[
i].block<4, 2>(0, 2);
8109 jcnt = jmax - jmin + 1;
8111 lp_bgn = forward ? nbremupstream : nbremupstream - 1;
8114 for (
int i = lp_bgn; forward ? (
i < lp_end) : (
i > lp_end);
i += (forward ? 1 : -1)) {
8116 i == bremno + (forward ? -1 : 1) &&
8118 prevstate->materialEffects() &&
8119 prevstate->materialEffects()->sigmaDeltaE() > 0
8121 jacbrem[
i].block(0, jmin, 4, jcnt) = jac.block(0, jmin, 4, jcnt);
8122 jacbrem[
i](4, 4) = jac(4, 4);
8124 calculateJac(jac, jacbrem[
i], jmin, jmax);
8128 Eigen::MatrixXd & derivmat = state->derivatives();
8129 int scatterPos = nperpars + 2 * nscats +
i;
8131 derivmat.block<5, 1>(0, scatterPos) = (forward ? .001 : -.001) * jacbrem[
i].block<5, 1>(0, 4);
8135 calculateJac(jac, jacvertex, 0, 4);
8139 Eigen::MatrixXd & derivmat = state->derivatives();
8140 derivmat.block(0, 0, 4, nperpars) = jacvertex.block(0, 0, 4, nperpars);
8142 if (nperpars == 5) {
8143 derivmat.col(4).segment(0, 4) *= .001;
8144 derivmat(4, 4) = .001 * jacvertex(4, 4);
8150 (!trajectory.prefit() ||
states[hitno]->materialEffects()->deltaE() == 0)
8152 scatno += (forward ? 1 : -1);
8156 states[hitno]->materialEffects() &&
8157 states[hitno]->materialEffects()->sigmaDeltaE() > 0
8159 bremno += (forward ? 1 : -1);
8162 prevstate =
states[hitno].get();
◆ calculateTrackErrors()
Definition at line 8169 of file GlobalChi2Fitter.cxx.
8177 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
8178 int nstatesupstream = trajectory.numberOfUpstreamStates();
8180 GXFTrackState *prevstate =
nullptr;
8181 int i = nstatesupstream;
8182 for (
int j = 0; j < (
int)
states.size(); j++) {
8183 if (j < nstatesupstream) {
8190 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
8191 if (stateno == 0 || stateno == nstatesupstream) {
8192 prevstate =
nullptr;
8195 std::unique_ptr<GXFTrackState> & state =
states[
index];
8196 if (state->materialEffects() !=
nullptr) {
8197 prevstate = state.get();
8201 if (!state->hasTrackCovariance()) {
8202 state->zeroTrackCovariance();
8204 AmgMatrix(5, 5) & trackerrmat = state->trackCovariance();
8206 if ((prevstate !=
nullptr) &&
8210 Eigen::Matrix<double, 5, 5> & jac = state->jacobian();
8213 trackerrmat = jac * prevcov * jac.transpose();
8217 trackerrmat = derivatives * fullcovmat * derivatives.transpose();
8221 const MeasurementBase *measurement = state->measurement();
8222 const Amg::MatrixX & meascov = measurement->localCovariance();
8224 ParamDefsAccessor paraccessor;
8228 bool errorok =
true;
8229 for (
int i = 0;
i < 5;
i++) {
8230 if (measurement->localParameters().contains(paraccessor.pardef[
i])) {
8232 && trackerrmat(
i,
i) > meascov(j, j)) {
8234 double scale = std::sqrt(meascov(j, j) / trackerrmat(
i,
i));
8235 trackerrmat(
i,
i) = meascov(j, j);
8236 for (
int k = 0;
k < 5;
k++) {
8246 for (
int i = 0;
i < 5;
i++) {
8250 for (
int j = 0; j < 5; j++) {
8257 if (trajectory.m_straightline) {
8258 trackerrmat(4, 4) = 1
e-20;
8262 state->trackParameters();
8264 std::optional<
AmgMatrix(5, 5)> trkerrmat;
8266 if (state->hasTrackCovariance()) {
8267 trkerrmat = (state->trackCovariance());
8269 trkerrmat = std::nullopt;
8272 const AmgVector(5) & tpars = tmptrackpar->parameters();
8273 std::unique_ptr<const TrackParameters> trackpar(
8274 tmptrackpar->associatedSurface().createUniqueTrackParameters(tpars[0],
8279 std::move(trkerrmat))
8281 state->setTrackParameters(std::move(trackpar));
8282 FitQualityOnSurface fitQual{};
8284 if (errorok && trajectory.nDOF() > 0) {
8285 fitQual =
m_updator->fullStateFitQuality(
8286 *state->trackParameters(),
8287 measurement->localParameters(),
8288 measurement->localCovariance()
8291 fitQual = FitQualityOnSurface(0, state->numberOfMeasuredParameters());
8294 state->setFitQuality(fitQual);
8296 prevstate = state.get();
◆ calculateTrackParameters()
Definition at line 7783 of file GlobalChi2Fitter.cxx.
7791 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
7792 int nstatesupstream = trajectory.numberOfUpstreamStates();
7793 const TrackParameters *prevtrackpar = trajectory.referenceParameters();
7794 std::unique_ptr<const TrackParameters> tmptrackpar;
7796 for (
int hitno = nstatesupstream - 1; hitno >= 0; hitno--) {
7800 DistanceSolution distsol = surf1.straightLineDistanceEstimate(
7801 prevtrackpar->position(), prevtrackpar->momentum().unit()
7808 distsol.numberOfSolutions() > 0 &&
7809 prevtrackpar != trajectory.referenceParameters()
7819 trajectory.m_fieldprop,
7826 (
rv.m_parameters !=
nullptr) &&
7827 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7833 if (
rv.m_parameters ==
nullptr) {
7834 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7835 " pos: " << prevtrackpar->position() <<
" destination surface: " << surf1);
7839 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7843 if (
rv.m_jacobian != std::nullopt) {
7845 states[hitno]->materialEffects() !=
nullptr &&
7846 states[hitno]->materialEffects()->deltaE() != 0 &&
7847 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7848 !trajectory.m_straightline
7850 double p = 1. / std::abs(currenttrackpar->parameters()[
Trk::qOverP]);
7851 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7852 double mass = trajectory.mass();
7853 double newp = std::sqrt(
p *
p + 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de);
7854 (*
rv.m_jacobian) (4, 4) = ((
p +
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7857 states[hitno]->setJacobian(*
rv.m_jacobian);
7858 }
else if (calcderiv) {
7863 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7865 if (meff !=
nullptr && hitno != 0) {
7866 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7867 surf, *meff, *
states[hitno]->trackParameters(), trajectory.mass(), -1
7870 if (std::holds_alternative<FitterStatusCode>(
r)) {
7871 return std::get<FitterStatusCode>(
r);
7874 tmptrackpar = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7875 prevtrackpar = tmptrackpar.get();
7877 prevtrackpar = currenttrackpar;
7881 prevtrackpar = trajectory.referenceParameters();
7883 for (
int hitno = nstatesupstream; hitno < (
int)
states.size(); hitno++) {
7886 DistanceSolution distsol = surf.straightLineDistanceEstimate(prevtrackpar->position(), prevtrackpar->momentum().unit());
7890 if (
distance < 0 && distsol.numberOfSolutions() > 0 && prevtrackpar != trajectory.referenceParameters()) {
7899 trajectory.m_fieldprop,
7905 (
rv.m_parameters !=
nullptr) &&
7907 (prevtrackpar->position() -
rv.m_parameters->position()).mag() > 5 *
mm
7912 if (
rv.m_parameters ==
nullptr) {
7913 ATH_MSG_DEBUG(
"propagation failed, prev par: " << *prevtrackpar <<
7914 " pos: " << prevtrackpar->
7915 position() <<
" destination surface: " << surf);
7919 if (
rv.m_jacobian != std::nullopt) {
7921 states[hitno]->materialEffects() !=
nullptr &&
7922 states[hitno]->materialEffects()->deltaE() != 0 &&
7923 states[hitno]->materialEffects()->sigmaDeltaE() <= 0 &&
7924 !trajectory.m_straightline
7926 double p = 1 / std::abs(
rv.m_parameters->parameters()[
Trk::qOverP]);
7927 double de = std::abs(
states[hitno]->materialEffects()->deltaE());
7928 double mass = trajectory.mass();
7929 double newp =
p *
p - 2 * de * std::sqrt(
mass *
mass +
p *
p) + de * de;
7932 newp = std::sqrt(newp);
7935 (*
rv.m_jacobian) (4, 4) = ((
p -
p * de / std::sqrt(
p *
p +
mass *
mass)) / newp) *
p *
p / (newp * newp);
7938 states[hitno]->setJacobian(*
rv.m_jacobian);
7939 }
else if (calcderiv) {
7944 GXFMaterialEffects *meff =
states[hitno]->materialEffects();
7946 if (meff !=
nullptr) {
7947 std::variant<std::unique_ptr<const TrackParameters>, FitterStatusCode>
r =
updateEnergyLoss(
7948 surf, *meff, *
rv.m_parameters, trajectory.mass(), +1
7951 if (std::holds_alternative<FitterStatusCode>(
r)) {
7952 return std::get<FitterStatusCode>(
r);
7955 rv.m_parameters = std::move(
std::get<std::unique_ptr<const TrackParameters>>(
r));
7958 states[hitno]->setTrackParameters(std::move(
rv.m_parameters));
7959 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 7757 of file GlobalChi2Fitter.cxx.
7766 PropagationResult
rv;
7769 ctx, prev,
ts, propdir, bf, calcderiv, holesearch
7772 if (
rv.m_parameters ==
nullptr) {
7773 propdir = invertPropdir(propdir);
7776 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 7718 of file GlobalChi2Fitter.cxx.
7727 std::unique_ptr<const TrackParameters>
rv;
7728 std::optional<TransportJacobian> jac{};
7739 if (
rv !=
nullptr && calcderiv) {
7744 std::optional<std::vector<std::unique_ptr<TrackParameters>>> extrapolation;
7750 return PropagationResult {
7753 std::move(extrapolation)
◆ compensatePhiWeights()
Definition at line 5959 of file GlobalChi2Fitter.cxx.
5964 const int nPerPars = trajectory.numberOfPerigeeParameters();
5965 std::size_t scatno = 0;
5967 for (
auto & state : trajectory.trackStates()) {
5968 const GXFMaterialEffects *meff = state->materialEffects();
5970 if (meff ==
nullptr || meff->sigmaDeltaPhi() == 0) {
5974 if (scatno >= cache.m_phiweight.size()) {
5976 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5977 throw std::range_error(
message.str());
5980 const bool isValidPlaneSurface =
5982 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5984 if (meff->deltaE() == 0 || isValidPlaneSurface) {
5985 const int scatNoIndex = 2 * scatno + nPerPars;
5986 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5987 cache.m_phiweight[scatno] = 1;
5994 if (meff->sigmaDeltaPhi() != 0) {
◆ correctAngles()
bool Trk::GlobalChi2Fitter::correctAngles |
( |
double & |
phi, |
|
|
double & |
theta |
|
) |
| |
|
staticprivate |
◆ ensureValidEntranceCalo()
bool Trk::GlobalChi2Fitter::ensureValidEntranceCalo |
( |
const EventContext & |
ctx, |
|
|
Cache & |
cache |
|
) |
| const |
|
private |
Definition at line 8531 of file GlobalChi2Fitter.cxx.
8532 if (cache.m_caloEntrance ==
nullptr) {
8536 cache.m_caloEntrance =
geometry->trackingVolume(
"InDet::Containers::InnerDetector");
8544 if (cache.m_caloEntrance ==
nullptr) {
8549 return cache.m_caloEntrance !=
nullptr;
◆ ensureValidEntranceMuonSpectrometer()
bool Trk::GlobalChi2Fitter::ensureValidEntranceMuonSpectrometer |
( |
const EventContext & |
ctx, |
|
|
Cache & |
cache |
|
) |
| const |
|
private |
Definition at line 8552 of file GlobalChi2Fitter.cxx.
8553 if (cache.m_msEntrance ==
nullptr) {
8557 cache.m_msEntrance =
geometry->trackingVolume(
"MuonSpectrometerEntrance");
8565 if (cache.m_msEntrance ==
nullptr) {
8570 return cache.m_msEntrance !=
nullptr;
◆ fillAfromMeasurements()
Definition at line 5756 of file GlobalChi2Fitter.cxx.
5761 const int nFitPars = trajectory.numberOfFitParameters();
5762 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5764 for (
int k = 0;
k < nFitPars;
k++) {
5765 for (
int l =
k;
l < nFitPars;
l++) {
5766 const int minMeas =
std::max(cache.m_firstmeasurement[
k], cache.m_firstmeasurement[
l]);
5767 const int maxMeas =
std::min(cache.m_lastmeasurement[
k], cache.m_lastmeasurement[
l]);
5770 for (
int measno = minMeas; measno < maxMeas; measno++) {
5771 a_kl += weightDeriv(measno,
k) * weightDeriv(measno,
l);
5774 a.fillSymmetric(
l,
k, a_kl);
◆ fillAfromScatterers()
Definition at line 5779 of file GlobalChi2Fitter.cxx.
5783 const int nFitPars = trajectory.numberOfFitParameters();
5784 const int nPerPars = trajectory.numberOfPerigeeParameters();
5785 const int nScatPars = 2 * trajectory.numberOfScatterers();
5786 const int nBrem = trajectory.numberOfBrems();
5787 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5790 const auto & scatSigmas = trajectory.scatteringSigmas();
5792 const int nMeas = (
int)
res.size();
5799 for (
int k = nPerPars;
k < nPerPars + nScatPars;
k += 2) {
5809 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5810 for (
int k = 4;
k < nFitPars;
k++) {
5812 k = nPerPars + nScatPars;
5815 for (
int l =
k;
l < nFitPars;
l++) {
5817 l = nPerPars + nScatPars;
5820 const double a_kl =
a(
l,
k) + weightDeriv(measno,
k) * weightDeriv(measno,
l);
5821 a.fillSymmetric(
l,
k, a_kl);
◆ fillBfromMeasurements()
Definition at line 5712 of file GlobalChi2Fitter.cxx.
5717 const int nFitPars = trajectory.numberOfFitParameters();
5718 const int nPerPars = trajectory.numberOfPerigeeParameters();
5719 const int nScatPars = 2 * trajectory.numberOfScatterers();
5720 const int nBrem = trajectory.numberOfBrems();
5721 const Amg::MatrixX & weightDeriv = trajectory.weightedResidualDerivatives();
5726 const int nMeas = (
int)
res.size();
5728 for (
int k = 0;
k < nFitPars;
k++) {
5729 const int minMeasK = cache.m_firstmeasurement[
k];
5730 const int maxMeasK = cache.m_lastmeasurement[
k];
5737 for (
int measno = minMeasK; measno < maxMeasK; measno++) {
5738 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
5748 if (
k == 4 ||
k >= nPerPars + nScatPars) {
5749 for (
int measno = nMeas - nBrem; measno < nMeas; measno++) {
5750 b[
k] +=
res[measno] * (1. /
error[measno]) * weightDeriv(measno,
k);
◆ fillDerivatives()
void Trk::GlobalChi2Fitter::fillDerivatives |
( |
GXFTrajectory & |
traj | ) |
const |
|
private |
Definition at line 5501 of file GlobalChi2Fitter.cxx.
5506 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5510 int nscatupstream = trajectory.numberOfUpstreamScatterers();
5511 int nbremupstream = trajectory.numberOfUpstreamBrems();
5512 int nscat = trajectory.numberOfScatterers();
5513 int nbrem = trajectory.numberOfBrems();
5514 int nperparams = trajectory.numberOfPerigeeParameters();
5516 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5519 int nmeas = (
int) weightderiv.rows();
5521 ParamDefsAccessor paraccessor;
5523 for (std::unique_ptr<GXFTrackState> & state :
states) {
5526 const MeasurementBase *measbase = state->measurement();
5527 const auto [scatmin, scatmax] = std::minmax(scatno, nscatupstream);
5528 const auto [bremmin, bremmax] = std::minmax(bremno, nbremupstream);
5535 const double sinStereo =
5537 state->sinStereo() :
5539 const double cosStereo =
5541 std::sqrt(1 -
std::pow(sinStereo, 2)) :
5549 auto getThisDeriv = [sinStereo, cosStereo, &derivatives](
int i,
int j) ->
double {
5550 if (
i == 0 && sinStereo != 0) {
5551 return derivatives(0, j) * cosStereo + sinStereo * derivatives(1, j);
5553 return derivatives(
i, j);
5557 for (
int i = 0;
i < 5;
i++) {
5558 if (!measbase->localParameters().contains(paraccessor.pardef[
i])) {
5569 if (trajectory.numberOfPerigeeParameters() > 0) {
5570 const int cols = trajectory.m_straightline ? 4 : 5;
5572 if (
i == 0 && sinStereo != 0) {
5573 weightderiv.row(measno).head(
cols) =
5574 (derivatives.row(0).head(
cols) * cosStereo +
5575 sinStereo * derivatives.row(1).head(
cols)) /
5578 weightderiv.row(measno).head(
cols) = derivatives.row(
i).head(
cols) /
error[measno];
5582 for (
int j = scatmin; j < scatmax; j++) {
5583 if (trajectory.prefit() == 1) {
5584 const int index = nperparams + j;
5587 const int index = nperparams + 2 * j;
5589 weightderiv(measno,
index + 1) = getThisDeriv(
i,
index + 1) /
error[measno];
5593 for (
int j = bremmin; j < bremmax; j++) {
5594 const int index = j + nperparams + 2 * nscat;
5601 double *
errors = state->measurementErrors();
5602 for (
int i = 0;
i < 5;
i++) {
5609 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5614 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5616 double qoverpbrem = limitInversePValue(1000 * state->trackParameters()->parameters()[
Trk::qOverP]);
5617 double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5619 double mass = .001 * trajectory.mass();
5621 const auto thisMeasurementIdx{nmeas - nbrem + bremno};
5625 auto multiplier = [] (
double mass,
double qOverP){
5628 const auto qoverpTerm {multiplier(
mass, qoverp) /
error[thisMeasurementIdx]};
5629 const auto qoverpBremTerm {multiplier(
mass, qoverpbrem) /
error[thisMeasurementIdx]};
5631 if (trajectory.numberOfPerigeeParameters() > 0) {
5632 weightderiv(thisMeasurementIdx, 4) = qoverpBremTerm - qoverpTerm;
5635 const auto bremNoBase = nperparams + 2 * nscat;
5636 if (bremno < nbremupstream) {
5637 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpTerm;
5638 for (
int bremno2 = bremno + 1; bremno2 < nbremupstream; bremno2++) {
5639 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpTerm - qoverpBremTerm;
5642 weightderiv(thisMeasurementIdx, bremNoBase + bremno) = qoverpBremTerm;
5643 for (
int bremno2 = nbremupstream; bremno2 < bremno; bremno2++) {
5644 weightderiv(thisMeasurementIdx, bremNoBase + bremno2) = qoverpBremTerm - qoverpTerm;
◆ fillFirstLastMeasurement()
void Trk::GlobalChi2Fitter::fillFirstLastMeasurement |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
trajectory |
|
) |
| const |
|
private |
Definition at line 5652 of file GlobalChi2Fitter.cxx.
5656 const int nFitPars = trajectory.numberOfFitParameters();
5657 const int nPerPars = trajectory.numberOfPerigeeParameters();
5658 const int nScatPars = 2 * trajectory.numberOfScatterers();
5659 const int nBrem = trajectory.numberOfBrems();
5660 const int nUpstreamStates = trajectory.numberOfUpstreamStates();
5663 const int nMeas = (
int)
res.size();
5665 cache.m_firstmeasurement.resize(nFitPars);
5666 cache.m_lastmeasurement.resize(nFitPars);
5668 for (
int i = 0;
i < nPerPars;
i++) {
5669 cache.m_firstmeasurement[
i] = 0;
5670 cache.m_lastmeasurement[
i] = nMeas - nBrem;
5676 for (
int i = 0;
i < (
int) trajectory.trackStates().size();
i++) {
5677 const std::unique_ptr<GXFTrackState> & state = trajectory.trackStates()[
i];
5678 const GXFMaterialEffects *meff = state->materialEffects();
5680 if (meff ==
nullptr) {
5681 measno += state->numberOfMeasuredParameters();
5685 const int firstMeasurement =
i < nUpstreamStates ? 0 : measno;
5686 const int lastMeasurement =
i < nUpstreamStates ? measno : nMeas - nBrem;
5688 if (meff->sigmaDeltaTheta() != 0
5689 && (trajectory.prefit() == 0 || meff->deltaE() == 0)) {
5690 const int scatterPos = nPerPars + 2 * scatno;
5692 cache.m_firstmeasurement[scatterPos] = firstMeasurement;
5693 cache.m_lastmeasurement[scatterPos] = lastMeasurement;
5695 cache.m_firstmeasurement[scatterPos + 1] = firstMeasurement;
5696 cache.m_lastmeasurement[scatterPos + 1] = lastMeasurement;
5701 if (meff->sigmaDeltaE() > 0) {
5702 const int bremPos = nPerPars + nScatPars + bremno;
5704 cache.m_firstmeasurement[bremPos] = firstMeasurement;
5705 cache.m_lastmeasurement[bremPos] = lastMeasurement;
◆ fillResidualsAndErrors()
Definition at line 5072 of file GlobalChi2Fitter.cxx.
5083 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
5103 const int nmeas = (
int)
res.size();
5104 const int nbrem = trajectory.numberOfBrems();
5105 const int nperpars = trajectory.numberOfPerigeeParameters();
5112 const int nidhits = trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits();
5113 const int nDOF = trajectory.nDOF();
5114 const bool doNewPseudoMeasurements = (
5118 std::abs((trajectory.prevchi2() - trajectory.chi2()) / nDOF) < 15 &&
5119 nidhits < trajectory.numberOfHits() &&
5120 (nperpars == 0 || nidhits > 0)
5131 double maxbrempull = -0.2;
5136 ParamDefsAccessor paraccessor;
5145 for (
int hitno = 0; hitno < (
int)
states.size(); hitno++) {
5146 std::unique_ptr<GXFTrackState> & state =
states[hitno];
5149 const MeasurementBase *measbase = state->measurement();
5160 doNewPseudoMeasurements &&
5162 !state->associatedSurface().isFree() &&
5163 !state->isRecalibrated()
5168 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
5171 currenttrackpar->associatedSurface()
5174 state->setMeasurement(std::move(newpseudo));
5175 measbase = state->measurement();
5182 double *
errors = state->measurementErrors();
5184 for (
int i = 0;
i < 5;
i++) {
5188 if (!measbase->localParameters().contains(paraccessor.pardef[
i])) {
5204 res[measno] = residuals[
i];
5219 double *
errors = state->measurementErrors();
5220 for (
int i = 0;
i < 5;
i++) {
5233 ((trajectory.prefit() == 0) || state->materialEffects()->deltaE() == 0)
5235 const double deltaPhi = state->materialEffects()->deltaPhi();
5236 const double measDeltaPhi = state->materialEffects()->measuredDeltaPhi();
5237 const double sigma2deltaPhi =
std::pow(state->materialEffects()->sigmaDeltaPhi(), 2);
5238 const double deltaTheta = state->materialEffects()->deltaTheta();
5239 const double sigma2deltaTheta =
std::pow(state->materialEffects()->sigmaDeltaTheta(), 2);
5241 if (trajectory.prefit() != 1) {
5242 b[nperpars + 2 * scatno] -= (
deltaPhi - measDeltaPhi) / sigma2deltaPhi;
5243 b[nperpars + 2 * scatno + 1] -= deltaTheta / sigma2deltaTheta;
5245 b[nperpars + scatno] -= deltaTheta / sigma2deltaTheta;
5250 deltaTheta * deltaTheta / sigma2deltaTheta
5259 if ((state->materialEffects() !=
nullptr) && state->materialEffects()->sigmaDeltaE() > 0) {
5260 double averagenergyloss = std::abs(state->materialEffects()->deltaE());
5262 const double qoverp = limitInversePValue(qoverpbrem - state->materialEffects()->delta_p());
5263 const double pbrem = 1. / std::abs(qoverpbrem);
5264 const double p = 1. / std::abs(qoverp);
5265 const double mass = .001 * trajectory.mass();
5267 const double bremEnergy = std::sqrt(pbrem * pbrem +
mass *
mass);
5269 const double resMaterial = .001 * averagenergyloss -
energy + bremEnergy;
5270 res[nmeas - nbrem + bremno] = resMaterial;
5272 const double sigde = state->materialEffects()->sigmaDeltaE();
5273 const double sigdepos = state->materialEffects()->sigmaDeltaEPos();
5274 const double sigdeneg = state->materialEffects()->sigmaDeltaENeg();
5276 double errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5277 error[nmeas - nbrem + bremno] = errorMaterial;
5288 if (state->materialEffects()->isKink()) {
5289 maxbrempull = -999999999;
5290 state_maxbrempull =
nullptr;
5294 cache.m_asymeloss &&
5296 trajectory.prefit() == 0 &&
5298 sigde != sigdepos &&
5301 const double elosspull = resMaterial / errorMaterial;
5303 if (trajectory.mass() > 100) {
5308 if (std::abs(elosspull) > 1) {
5309 if (elosspull < -1) {
5310 state->materialEffects()->setSigmaDeltaE(sigdepos);
5312 state->materialEffects()->setSigmaDeltaE(sigdeneg);
5315 errorMaterial = .001 * state->materialEffects()->sigmaDeltaE();
5316 error[nmeas - nbrem + bremno] = errorMaterial;
5318 }
else if ((trajectory.numberOfTRTHits() == 0) ||
it >= 3) {
5328 !state->materialEffects()->isKink() && (
5329 (
m_fixbrem == -1 && elosspull < maxbrempull) ||
5333 bremno_maxbrempull = bremno;
5334 state_maxbrempull = state.get();
5335 maxbrempull = elosspull;
5344 trajectory.prefit() == 0 &&
5345 state->materialEffects()->sigmaDeltaPhi() == 0 &&
5346 state->materialEffects()->isMeasuredEloss() &&
5347 resMaterial / (.001 * state->materialEffects()->sigmaDeltaEAve()) > 2.5
5352 std::vector<MaterialEffectsOnTrack> calomeots =
5357 parforcalo->associatedSurface(),
5365 if (calomeots.size() == 3) {
5366 averagenergyloss = std::abs(calomeots[1].energyLoss()->deltaE());
5367 const double newres = .001 * averagenergyloss -
energy + bremEnergy;
5368 const double newerr = .001 * calomeots[1].energyLoss()->sigmaDeltaE();
5370 const double oldPull = resMaterial / errorMaterial;
5371 const double newPull = newres / newerr;
5373 if (std::abs(newPull) < std::abs(oldPull)) {
5374 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
5376 state->materialEffects()->setEloss(std::unique_ptr<EnergyLoss>(calomeots[1].energyLoss()->
clone()));
5377 state->materialEffects()->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
5378 res[nmeas - nbrem + bremno] = newres;
5379 error[nmeas - nbrem + bremno] = newerr;
5383 state->materialEffects()->setMeasuredEloss(
false);
5393 for (
int imeas = 0; imeas < nmeas; imeas++) {
5394 if (
error[imeas] == 0) {
5404 trajectory.setPrevChi2(trajectory.chi2());
5405 trajectory.setChi2(
chi2);
◆ finalize()
StatusCode Trk::GlobalChi2Fitter::finalize |
( |
| ) |
|
|
overridevirtual |
Definition at line 290 of file GlobalChi2Fitter.cxx.
293 if (m_fit_status[
S_FITS] > 0) {
296 <<
" track fits failed because of a matrix inversion failure");
298 <<
" tracks were rejected by the outlier logic");
300 <<
" track fits failed because of a propagation failure");
302 <<
" track fits failed because of an invalid angle (theta/phi)");
304 <<
" track fits failed because the fit did not converge");
306 <<
" tracks did not pass the chi^2 cut");
308 <<
" tracks were killed by the energy loss update");
311 return StatusCode::SUCCESS;
◆ fit() [1/6]
Definition at line 2401 of file GlobalChi2Fitter.cxx.
2408 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Meas'BaseSet,,)");
2413 GXFTrajectory trajectory;
2416 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2421 for (
const auto *itSet : rots) {
2422 if (itSet ==
nullptr) {
2423 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2429 std::unique_ptr<const TrackParameters> startpar(param.clone());
2432 matEffects ==
muon &&
2433 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0
2435 cache.m_matfilled =
true;
2436 trajectory.setPrefit(2);
2438 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2440 cache.m_matfilled =
false;
2442 if (!trajectory.converged()) {
2446 trajectory.setConverged(
false);
2447 const TrackParameters *firstpar = trajectory.trackStates()[0]->trackParameters();
2448 const TrackParameters *lastpar = trajectory.trackStates().back()->trackParameters();
2450 PerigeeSurface persurf(firstpar->position() - 10 * firstpar->momentum().unit());
2456 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2459 firstpar->associatedSurface()
2462 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2469 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2472 lastpar->associatedSurface()
2475 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2478 if (!trajectory.m_straightline) {
2479 trajectory.setPrefit(3);
2480 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2481 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2482 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2487 myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects);
2489 cache.m_matfilled =
true;
2491 if (!trajectory.converged()) {
2496 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2497 startpar = trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2498 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2502 trajectory.setPrefit(0);
2505 firstpar = trajectory.trackStates().front()->trackParameters();
2510 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2513 firstpar->associatedSurface()
2516 trajectory.trackStates().front()->setMeasurement(std::move(newpseudo));
2520 trajectory.trackStates().front()->setMeasurementErrors(
errors);
2524 lastpar = trajectory.trackStates().back()->trackParameters();
2529 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2532 lastpar->associatedSurface()
2535 trajectory.trackStates().back()->setMeasurement(std::move(newpseudo));
2539 trajectory.trackStates().back()->setMeasurementErrors(
errors);
2543 std::unique_ptr<Track>
track;
2545 if (startpar !=
nullptr) {
2546 track.reset(
myfit(ctx,cache, trajectory, *startpar, runOutlier, matEffects));
2549 if (track !=
nullptr) {
2552 cache.m_matfilled =
false;
◆ fit() [2/6]
Definition at line 2180 of file GlobalChi2Fitter.cxx.
2187 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(PRDS,TP,)");
2190 for (
const auto *prd : prds) {
2191 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2193 const PlaneSurface *plsurf =
nullptr;
2196 plsurf =
static_cast < const PlaneSurface *
>(&prdsurf);
2198 const StraightLineSurface *slsurf =
nullptr;
2201 slsurf =
static_cast < const StraightLineSurface *
>(&prdsurf);
2203 if ((slsurf ==
nullptr) && (plsurf ==
nullptr)) {
2204 ATH_MSG_ERROR(
"Surface is neither PlaneSurface nor StraightLineSurface!");
2209 }
else if (slsurf !=
nullptr) {
2218 }
else if (plsurf !=
nullptr) {
2219 if (param.covariance() !=
nullptr) {
2241 if (rot !=
nullptr) {
2242 rots.push_back(rot);
2246 std::unique_ptr<Track>
track =
2247 fit(ctx, rots, param, runOutlier, matEffects);
2249 for (
const auto *rot : rots) {
◆ fit() [3/6]
Definition at line 2257 of file GlobalChi2Fitter.cxx.
2264 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Meas'BaseSet,,)");
2269 GXFTrajectory trajectory;
2272 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
2279 if (minpar ==
nullptr) {
2280 minpar = *(inputTrack.trackParameters()->begin());
2289 bool old_reintoutl = cache.m_reintoutl;
2290 cache.m_reintoutl =
false;
2291 bool tmpasymeloss = cache.m_asymeloss;
2294 cache.m_asymeloss =
true;
2297 for (; itStates != endState; ++itStates) {
2301 !trajectory.trackStates().empty() &&
2302 (trajectory.trackStates().back()->materialEffects() !=
nullptr) &&
2303 trajectory.trackStates().back()->materialEffects()->sigmaDeltaE() > 50.001
2305 trajectory.trackStates().back()->materialEffects()->setKink(
true);
2309 cache.m_reintoutl = old_reintoutl;
2311 for (
const auto & measBase : addMeasColl) {
2312 if (measBase ==
nullptr) {
2313 ATH_MSG_WARNING(
"There is an empty MeasurementBase object in the track! Skip this object..");
2321 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2322 cache.m_asymeloss = tmpasymeloss;
2324 if (track !=
nullptr) {
2326 inputTrack.fitQuality()->numberDoF() != 0 ?
2327 inputTrack.fitQuality()->chiSquared() / inputTrack.fitQuality()->numberDoF() :
2331 track->fitQuality()->numberDoF() != 0 ?
2332 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() :
2335 if (
m_extensioncuts && runOutlier && newqual > 2 && newqual > 2 * oldqual) {
2339 track.reset(
nullptr);
2343 if (track !=
nullptr) {
◆ fit() [4/6]
Definition at line 2353 of file GlobalChi2Fitter.cxx.
2360 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,PRDS,)");
2364 for (
const auto *prd : prds) {
2365 const Surface & prdsurf = (*prd).detectorElement()->surface((*prd).identify());
2368 std::unique_ptr<const TrackParameters>trackparForCorrect(
2369 hitparam->associatedSurface().createUniqueTrackParameters(
2384 rot =
m_ROTcreator->correct(*prd, *trackparForCorrect, ctx);
2387 if (rot !=
nullptr) {
2388 rots.push_back(rot);
2392 std::unique_ptr<Track>
track =
fit(ctx,intrk, rots, runOutlier, matEffects);
2394 for (
const auto *rot : rots) {
◆ fit() [5/6]
Definition at line 1795 of file GlobalChi2Fitter.cxx.
1801 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,)");
1806 GXFTrajectory trajectory;
1809 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1814 return std::unique_ptr<Track>(
1815 fitIm(ctx, cache, inputTrack, runOutlier, matEffects));
◆ fit() [6/6]
Definition at line 317 of file GlobalChi2Fitter.cxx.
324 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,Track,)");
329 GXFTrajectory trajectory;
331 trajectory.m_straightline = (
332 !cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn()
339 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
340 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
342 bool measphi =
false;
344 for (
const auto *
i : *(muontrack->measurementsOnTrack())) {
349 rot = &crot->rioOnTrack(0);
356 const Surface *surf = &rot->associatedSurface();
359 double dotprod2 = measdir.dot(
361 surf->center().perp());
362 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
370 auto [firstidpar, lastidpar] = getFirstLastIdPar(*indettrack);
372 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
376 std::unique_ptr<const TrackParameters> parforcalo =
unique_clone(firstismuon ? firstidpar : lastidpar);
378 if (!cache.m_field_cache.solenoidOn()) {
379 const AmgVector(5) & newpars = parforcalo->parameters();
381 parforcalo=parforcalo->associatedSurface().createUniqueTrackParameters(
382 newpars[0], newpars[1], newpars[2], newpars[3], 1 / 5000., std::nullopt
386 std::vector < MaterialEffectsOnTrack > calomeots;
389 calomeots =
m_calotool->extrapolationSurfacesAndEffects(
393 parforcalo->associatedSurface(),
406 if (calomeots.empty()) {
411 std::unique_ptr<Track>
track;
414 cache.m_calomat =
false;
415 bool tmp2 = cache.m_extmat;
416 bool tmp4 = cache.m_idmat;
421 double qoverpid = measperid !=
nullptr ? measperid->parameters()[
Trk::qOverP] : 0;
422 double qoverpmuon = measpermuon !=
nullptr ? measpermuon->parameters()[
Trk::qOverP] : 0;
424 const AmgSymMatrix(5) * errmatid = measperid !=
nullptr ? measperid->covariance() :
nullptr;
425 const AmgSymMatrix(5) * errmatmuon = measpermuon !=
nullptr ? measpermuon->covariance() :
nullptr;
429 (errmatid !=
nullptr) &&
430 (errmatmuon !=
nullptr) &&
435 double piderror = std::sqrt((*errmatid) (4, 4)) / (qoverpid * qoverpid);
436 double pmuonerror = std::sqrt((*errmatmuon) (4, 4)) / (qoverpmuon * qoverpmuon);
437 double energyerror = std::sqrt(
438 calomeots[1].energyLoss()->sigmaDeltaE() *
439 calomeots[1].energyLoss()->sigmaDeltaE() + piderror * piderror +
440 pmuonerror * pmuonerror
444 (std::abs(calomeots[1].energyLoss()->deltaE()) -
445 std::abs(1 / qoverpid) + std::abs(1 / qoverpmuon)
448 ATH_MSG_DEBUG(
"Changing from measured to parametrized energy loss");
453 parforcalo->associatedSurface(),
458 if (calomeots.empty()) {
465 int nfits = cache.m_fit_status[
S_FITS];
466 bool firstfitwasattempted =
false;
469 if (!caloEntranceIsValid) {
474 (!cache.m_field_cache.toroidOn() && !cache.m_field_cache.solenoidOn()) ||
476 cache.m_getmaterialfromtrack &&
480 qoverpid * qoverpmuon > 0
485 if (cache.m_fit_status[
S_FITS] == (
unsigned int) (nfits + 1)) {
486 firstfitwasattempted =
true;
491 (track ==
nullptr) &&
492 !firstfitwasattempted &&
493 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn())
496 GXFTrajectory trajectory2;
497 trajectory2.m_straightline = trajectory.m_straightline;
498 trajectory2.m_fieldprop = trajectory.m_fieldprop;
499 trajectory = trajectory2;
503 bool pseudoupdated =
false;
505 if (track !=
nullptr) {
506 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
507 if (pseudostate ==
nullptr) {
518 if ((pseudostate ==
nullptr) || pseudostate->fitQuality().chiSquared() < 10) {
523 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
525 pseudostate->measurement()->localParameters(),
526 pseudostate->measurement()->localCovariance()
529 if (updpar ==
nullptr) {
536 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
541 pseudopar->associatedSurface()
544 pseudostate->setMeasurement(std::move(newpseudo));
548 pseudostate->setMeasurementErrors(
errors);
549 pseudoupdated =
true;
553 trajectory.setConverged(
false);
554 cache.m_matfilled =
true;
560 *
track->perigeeParameters(),
562 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
565 cache.m_matfilled =
false;
569 cache.m_fit_status[
S_FITS] = nfits + 1;
571 if (track !=
nullptr) {
572 track->info().addPatternReco(intrk1.info());
573 track->info().addPatternReco(intrk2.info());
577 cache.m_calomat =
tmp;
578 cache.m_extmat =
tmp2;
579 cache.m_idmat = tmp4;
◆ fitIm()
Definition at line 1849 of file GlobalChi2Fitter.cxx.
1857 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::fit(Track,,)");
1859 GXFTrajectory trajectory;
1862 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
1867 if (inputTrack.trackStateOnSurfaces()->empty()) {
1872 if (inputTrack.trackParameters()->empty()) {
1873 ATH_MSG_WARNING(
"Track without track parameters, cannot perform fit");
1877 std::unique_ptr<const TrackParameters> minpar =
unique_clone(inputTrack.perigeeParameters());
1881 if (minpar ==
nullptr) {
1882 minpar =
unique_clone(*(inputTrack.trackParameters()->begin()));
1885 bool tmpgetmat = cache.m_getmaterialfromtrack;
1891 cache.m_getmaterialfromtrack =
false;
1897 trajectory.trackStates().reserve(inputTrack.trackStateOnSurfaces()->size());
1899 const Surface *firsthitsurf =
nullptr;
1900 const Surface *lasthitsurf =
nullptr;
1902 bool hasmuon =
false;
1903 bool iscombined =
false;
1904 bool seenphimeas =
false;
1908 for (; itStates != endState; ++itStates) {
1909 const auto *
const pMeasurement = (**itStates).measurementOnTrack();
1911 (pMeasurement ==
nullptr) &&
1912 ((**itStates).materialEffectsOnTrack() !=
nullptr) &&
1913 matEffects != inputTrack.info().particleHypothesis()
1918 if (pMeasurement !=
nullptr) {
1919 const Surface *surf = &pMeasurement->associatedSurface();
1920 Identifier hitid = surf->associatedDetectorElementIdentifier();
1924 hitid = crot->rioOnTrack(0).identify();
1928 if (firsthitsurf ==
nullptr) {
1929 firsthitsurf = surf;
1937 if ((**itStates).trackParameters() !=
nullptr) {
1938 lastidpar = (**itStates).trackParameters();
1939 if (firstidpar ==
nullptr) {
1940 firstidpar = lastidpar;
1945 Amg::Vector3D measdir = surf->transform().rotation().col(0);
1947 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
1948 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
1950 if (std::abs(surf->center().z()) > 13000) {
1953 if (surf->center().perp() > 9000 && std::abs(surf->center().z()) < 13000) {
1965 if (iscombined && seenphimeas && (phiem || phibo)) {
1975 cache.m_getmaterialfromtrack && trajectory.numberOfScatterers() != 0 &&
1976 (hasmuon || cache.m_acceleration)
1978 cache.m_matfilled =
true;
1981 if (firstidpar == lastidpar) {
1982 firstidpar = lastidpar =
nullptr;
1987 !cache.m_matfilled &&
1989 m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier()) !=
1992 (firstidpar !=
nullptr)
1994 if (
m_DetID->
is_indet(firsthitsurf->associatedDetectorElementIdentifier())) {
2001 bool tmpacc = cache.m_acceleration;
2003 bool tmpsirecal = cache.m_sirecal;
2004 std::unique_ptr<Track> tmptrack =
nullptr;
2008 cache.m_fiteloss =
true;
2009 cache.m_sirecal =
false;
2012 cache.m_asymeloss =
true;
2015 tmptrack.reset(
myfit(ctx, cache, trajectory, *minpar,
false, matEffects));
2016 cache.m_sirecal = tmpsirecal;
2018 if (tmptrack ==
nullptr) {
2019 cache.m_matfilled =
false;
2020 cache.m_getmaterialfromtrack = tmpgetmat;
2021 cache.m_acceleration = tmpacc;
2022 cache.m_fiteloss = tmpfiteloss;
2027 bool isbrem =
false;
2029 unsigned int n_brem=0;
2031 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
2032 GXFMaterialEffects *meff = state->materialEffects();
2034 if (meff !=
nullptr) {
2038 const MaterialProperties *matprop = meff->materialProperties();
2040 double p = 1. / std::abs(layerpars->parameters()[
Trk::qOverP] - .0005 * meff->delta_p());
2042 std::optional<Amg::Vector2D> locpos(state->associatedSurface().globalToLocal(layerpars->position()));
2043 const Amg::Vector3D layerNormal(state->associatedSurface().normal(*locpos));
2044 double costracksurf = 1.;
2046 costracksurf = std::abs(layerNormal.dot(layerpars->momentum().unit()));
2048 double oldde = meff->deltaE();
2050 std::unique_ptr<EnergyLoss> eloss;
2051 double sigmascat = 0;
2053 if (matprop !=
nullptr) {
2054 eloss = std::make_unique<EnergyLoss>(
2055 m_elosstool->energyLoss(*matprop,
p, 1. / costracksurf,
2057 sigmascat = std::sqrt(
m_scattool->sigmaSquare(*matprop,
p, 1. / costracksurf, matEffects));
2059 if (eloss !=
nullptr) {
2060 meff->setDeltaE(eloss->deltaE());
2063 MaterialProperties tmpprop(1., meff->x0(), 0., 0., 0., 0.);
2064 sigmascat = std::sqrt(
m_scattool->sigmaSquare(tmpprop,
p, 1. / costracksurf, matEffects));
2067 meff->setScatteringSigmas(
2075 meff->setDeltaE(oldde);
2076 if (!meff->isKink()) {
2077 meff->setSigmaDeltaE(0);
2080 bremdp = meff->delta_p();
2082 }
else if (eloss !=
nullptr) {
2083 meff->setSigmaDeltaE(eloss->sigmaDeltaE());
2085 if (meff->sigmaDeltaE() > 0) {
2091 const AmgVector(5) & refpars = trajectory.referenceParameters()->parameters();
2092 minpar=trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
2093 refpars[0], refpars[1], refpars[2], refpars[3], refpars[4], std::nullopt
2097 cache.m_matfilled =
true;
2101 trajectory.brems().clear();
2103 trajectory.brems().resize(1);
2104 trajectory.brems()[0] = bremdp;
2107 cache.m_asymeloss =
false;
2108 trajectory.setNumberOfScatterers(nscats);
2111 trajectory.setNumberOfBrems(n_brem);
2115 std::unique_ptr<Track>
track(
myfit(ctx, cache, trajectory, *minpar, runOutlier, matEffects));
2117 bool pseudoupdated =
false;
2119 if ((track !=
nullptr) && hasid && hasmuon) {
2120 for (std::unique_ptr<GXFTrackState> & pseudostate : trajectory.trackStates()) {
2122 (pseudostate ==
nullptr) ||
2124 pseudostate->fitQuality().chiSquared() < 10
2130 std::unique_ptr<const TrackParameters> updpar(
m_updator->removeFromState(
2132 pseudostate->measurement()->localParameters(),
2133 pseudostate->measurement()->localCovariance()
2136 if (updpar ==
nullptr) {
2143 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
2146 pseudopar->associatedSurface()
2149 pseudostate->setMeasurement(std::move(newpseudo));
2153 pseudostate->setMeasurementErrors(
errors);
2154 pseudoupdated =
true;
2157 if (pseudoupdated) {
2158 trajectory.setConverged(
false);
2159 cache.m_matfilled =
true;
2161 cache.m_matfilled =
false;
2165 cache.m_matfilled =
false;
2166 cache.m_getmaterialfromtrack = tmpgetmat;
2167 cache.m_acceleration = tmpacc;
2168 cache.m_fiteloss = tmpfiteloss;
2170 if (track !=
nullptr) {
2172 const TrackInfo& old_info = inputTrack.info();
2173 track->info().addPatternReco(old_info);
2176 return track.release();
◆ holesearchExtrapolation()
Helper method which performs an extrapolation with additional logic for hole search.
This method is a wrapper around extrapolateStepwise from the extrapolator interface, with the added functionality that it will null any returned track parameters which are on the start and end surface.
- Parameters
-
[in] | ctx | An event context for extrapolation. |
[in] | src | The track parameters to start extrapolating from. |
[in] | dst | The track state to extrapolate to. |
[in] | propdir | The propagation direction. |
- Returns
- A vector of track states, just like normal extrapolation.
Definition at line 7667 of file GlobalChi2Fitter.cxx.
7677 std::vector<std::unique_ptr<TrackParameters>>
rv =
m_extrapolator->extrapolateStepwise(
7678 ctx,
src, dst.associatedSurface(), propdir,
false
7690 &
rv.front()->associatedSurface() == &dst.associatedSurface() ||
7691 &
rv.front()->associatedSurface() == &
src.associatedSurface() ||
7692 trackParametersClose(*
rv.front(),
src, 0.001) ||
7693 trackParametersClose(*
rv.front(), *dst.trackParameters(), 0.001)
7696 rv.front().reset(
nullptr);
7706 &
rv.back()->associatedSurface() == &dst.associatedSurface() ||
7707 &
rv.back()->associatedSurface() == &
src.associatedSurface() ||
7708 trackParametersClose(*
rv.back(),
src, 0.001) ||
7709 trackParametersClose(*
rv.back(), *dst.trackParameters(), 0.001)
7712 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 7210 of file GlobalChi2Fitter.cxx.
7223 for (
const std::unique_ptr<TrackParameters> &
tp : hc) {
7229 if (
tp ==
nullptr) {
7238 const TrkDetElementBase * de =
tp->associatedSurface().associatedDetectorElement();
7240 if (de ==
nullptr) {
7251 if (id_set.find(
id) != id_set.end()) {
7300 if (sct_set.find(
os) != sct_set.end()) {
7301 ++
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 7395 of file GlobalChi2Fitter.cxx.
7409 constexpr
uint min_meas = 3;
7410 if (std::count_if(
states.begin(),
states.end(), [](
const GXFTrackState &
s){ return s.getStateType(TrackStateOnSurface::Measurement); }) < min_meas) {
7414 bool seen_meas =
false;
7416 std::set<Identifier> id_set;
7417 std::set<Identifier> sct_set;
7423 for (std::size_t
i = 0;
i <
states.size() - 1;
i++) {
7446 double dist = (
beg.trackParameters()->position() -
end.trackParameters()->position()).
norm();
7448 bool zStartValid = std::abs(
beg.trackParameters()->position().z())<10000.;
7450 ATH_MSG_DEBUG(
"Pathological track parameter well outside of detector");
7451 ATH_MSG_DEBUG(
"Propagator might have issue with this, skipping");
7461 if (seen_meas && dist >= 2.5 && zStartValid) {
7468 std::optional<std::vector<std::unique_ptr<TrackParameters>>> & hc =
beg.getHoles();
7469 std::vector<std::unique_ptr<TrackParameters>>
states;
7477 if (hc.has_value()) {
7498 GXFTrackState & last =
states.back();
7509 last.trackParameters() !=
nullptr &&
7516 std::vector<std::unique_ptr<Trk::TrackParameters>> bl =
m_extrapolator->extrapolateBlindly(
7518 *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 7320 of file GlobalChi2Fitter.cxx.
7328 GXFTrackState * lastmeas =
nullptr;
7330 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7342 std::vector<std::reference_wrapper<GXFTrackState>>
rv;
7343 rv.reserve(trajectory.trackStates().size());
7350 for (
const std::unique_ptr<GXFTrackState> &
s : trajectory.trackStates()) {
7364 rv.emplace_back(*
s);
7371 const TrkDetElementBase * de =
s->trackParameters()->associatedSurface().associatedDetectorElement();
7373 if (de !=
nullptr) {
7386 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 8507 of file GlobalChi2Fitter.cxx.
8517 if (cond_obj ==
nullptr) {
8522 cond_obj->getInitializedCache(cache.m_field_cache);
◆ initialize()
StatusCode Trk::GlobalChi2Fitter::initialize |
( |
| ) |
|
|
overridevirtual |
Definition at line 211 of file GlobalChi2Fitter.cxx.
231 ATH_MSG_ERROR(
"Hole search requested but no boundary check tool provided.");
232 return StatusCode::FAILURE;
257 ATH_MSG_WARNING(
"FillDerivativeMatrix option selected, switching off acceleration!");
281 ATH_MSG_ERROR(
"Hole search requested but track summaries are disabled.");
282 return StatusCode::FAILURE;
287 return StatusCode::SUCCESS;
◆ isMuonTrack()
bool Trk::GlobalChi2Fitter::isMuonTrack |
( |
const Track & |
intrk1 | ) |
const |
|
private |
Definition at line 8465 of file GlobalChi2Fitter.cxx.
8466 const auto *pDataVector = intrk1.measurementsOnTrack();
8467 auto nmeas1 = pDataVector->size();
8468 const auto *pLastValue = (*pDataVector)[nmeas1 - 1];
8476 testrot =
static_cast<const RIO_OnTrack *
>(pLastValue);
8478 if (lastMeasIsCompetingRIO){
8480 testrot = &testcrot->rioOnTrack(0);
8484 if (testrot ==
nullptr) {
8485 const auto *pPenultimate = (*pDataVector)[nmeas1 - 2];
8488 if(penultimateIsRIO){
8489 testrot =
static_cast<const RIO_OnTrack *
>(pPenultimate);
8491 if (penultimateIsCompetingRIO){
8493 testrot = &testcrot->rioOnTrack(0);
8500 (testrot !=
nullptr) &&
◆ iterationsOfLastFit()
int Trk::GlobalChi2Fitter::iterationsOfLastFit |
( |
| ) |
const |
|
privatevirtual |
◆ mainCombinationStrategy()
Definition at line 583 of file GlobalChi2Fitter.cxx.
591 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::mainCombinationStrategy");
596 const Track *indettrack = firstismuon ? &intrk2 : &intrk1;
597 const Track *muontrack = firstismuon ? &intrk1 : &intrk2;
599 auto [tmpfirstidpar, tmplastidpar] = getFirstLastIdPar(*indettrack);
600 std::unique_ptr<const TrackParameters> firstidpar =
unique_clone(tmpfirstidpar);
601 std::unique_ptr<const TrackParameters> lastidpar =
unique_clone(tmplastidpar);
603 if ((firstidpar ==
nullptr) || (lastidpar ==
nullptr)) {
607 if (muontrack->trackStateOnSurfaces()->empty()) {
613 muontrack->trackStateOnSurfaces()->end() - 1 :
614 muontrack->trackStateOnSurfaces()->begin();
616 const MeasurementBase *closestmuonmeas =
nullptr;
617 std::unique_ptr<const TrackParameters> tp_closestmuon =
nullptr;
619 while (closestmuonmeas ==
nullptr) {
620 closestmuonmeas =
nullptr;
623 if ((**tsosit).measurementOnTrack() !=
nullptr) {
624 closestmuonmeas = (**tsosit).measurementOnTrack();
626 if (thispar !=
nullptr) {
627 const AmgVector(5) & parvec = thispar->parameters();
628 tp_closestmuon=thispar->associatedSurface().createUniqueTrackParameters(
629 parvec[0], parvec[1], parvec[2], parvec[3], parvec[4], std::nullopt
643 std::unique_ptr<const TrackParameters> tmppar;
646 if ((tp_closestmuon !=
nullptr) && msEntranceIsValid) {
648 ctx, *tp_closestmuon, *cache.m_msEntrance, propdir,
nonInteracting);
651 std::unique_ptr<const std::vector<const TrackStateOnSurface *>> matvec;
653 if (tmppar !=
nullptr) {
654 const Surface & associatedSurface = tmppar->associatedSurface();
655 std::unique_ptr<Surface> muonsurf =
nullptr;
659 const CylinderBounds *cylbounds =
static_cast <const CylinderBounds *
>(&associatedSurface.bounds());
661 double radius = cylbounds->r();
662 double hlength = cylbounds->halflengthZ();
663 muonsurf = std::make_unique<CylinderSurface>(trans,
radius + 1, hlength);
668 associatedSurface.center().z() > 0 ?
669 associatedSurface.center().z() + 1 :
670 associatedSurface.center().z() - 1
674 associatedSurface.center().x(),
675 associatedSurface.center().y(),
679 trans.translation() << newpos;
681 const DiscBounds *discbounds =
static_cast<const DiscBounds *
>(&associatedSurface.bounds());
682 double rmin = discbounds->rMin();
683 double rmax = discbounds->rMax();
684 muonsurf = std::make_unique<DiscSurface>(trans, rmin, rmax);
688 if (muonsurf !=
nullptr) {
700 std::vector<const TrackStateOnSurface *> tmp_matvec;
702 if ((matvec !=
nullptr) && !matvec->empty()) {
703 tmp_matvec = *matvec;
704 delete tmp_matvec.back();
705 tmp_matvec.pop_back();
707 for (
auto &
i : tmp_matvec) {
712 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(
i->materialEffectsOnTrack());
714 const Surface *matsurf = &meff->associatedSurface();
721 trajectory.m_fieldprop,
726 if (tmppar ==
nullptr) {
734 trajectory.m_fieldprop,
740 if (tmppar ==
nullptr) {
748 double de = std::abs(meff->energyLoss()->deltaE());
750 double newp2 = oldp * oldp + (!firstismuon ? 2 : -2) * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
757 tp_closestmuon=tmppar->associatedSurface().createUniqueTrackParameters(
758 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
774 for (; itStates != endState; ++itStates) {
779 bool tmpgetmat = cache.m_getmaterialfromtrack;
781 if ((*itStates)->materialEffectsOnTrack() !=
nullptr) {
783 cache.m_extmat =
false;
785 cache.m_idmat =
false;
788 const auto *
const pBaseMEOT = (*itStates)->materialEffectsOnTrack();
792 const auto *
const pMEOT =
static_cast<const MaterialEffectsOnTrack *
>((*itStates)->materialEffectsOnTrack());
793 if ((pMEOT->scatteringAngles() ==
nullptr) or (pMEOT->energyLoss() ==
nullptr)) {
794 cache.m_getmaterialfromtrack =
true;
800 cache.m_getmaterialfromtrack = tmpgetmat;
807 trajectory.trackStates().back()->setTrackParameters(
nullptr);
810 std::unique_ptr<const TrackParameters> firstscatpar;
811 std::unique_ptr<const TrackParameters> lastscatpar;
814 double newqoverpid = 0;
817 double de = std::abs(calomeots[1].energyLoss()->deltaE());
818 double sigmade = std::abs(calomeots[1].energyLoss()->sigmaDeltaE());
820 double pbefore = std::abs(1 / firstidpar->parameters()[
Trk::qOverP]);
821 double pafter = std::abs(1 / tp_closestmuon->parameters()[
Trk::qOverP]);
822 double elosspull = (pbefore - pafter - de) / sigmade;
824 if (std::abs(elosspull) > 10) {
825 if (elosspull > 10) {
826 newqoverpid = 1 / (de + pafter + 10 * sigmade);
828 newqoverpid = 1 / (de + pafter - 10 * sigmade);
831 if (tp_closestmuon->parameters()[
Trk::qOverP] * newqoverpid < 0) {
835 const AmgVector(5) & newpar = firstidpar->parameters();
836 firstidpar=firstidpar->associatedSurface().createUniqueTrackParameters(
837 newpar[0], newpar[1], newpar[2], newpar[3], newqoverpid, std::nullopt
845 if (lastidpar ==
nullptr) {
851 *(firstismuon ? tp_closestmuon.get() : lastidpar.get()),
852 calomeots[0].associatedSurface(),
855 trajectory.m_fieldprop,
859 if (firstscatpar ==
nullptr) {
865 *(firstismuon ? firstidpar : tp_closestmuon),
866 calomeots[2].associatedSurface(),
869 trajectory.m_fieldprop,
873 if (lastscatpar ==
nullptr) {
877 std::optional<TransportJacobian> jac1, jac2;
878 std::unique_ptr<const TrackParameters> elosspar;
880 double firstscatphi = 0;
881 double secondscatphi = 0;
882 double firstscattheta = 0;
883 double secondscattheta = 0;
884 double muonscatphi = 0;
885 double muonscattheta = 0;
887 const TrackParameters *idscatpar = !firstismuon ? firstscatpar.get() : lastscatpar.get();
888 const TrackParameters *muonscatpar = firstismuon ? firstscatpar.get() : lastscatpar.get();
890 newqoverpid = idscatpar->parameters()[
Trk::qOverP];
892 Amg::Vector3D calosegment = lastscatpar->position() - firstscatpar->position();
895 muonscattheta = calosegment.theta() - muonscatpar->parameters()[
Trk::theta];
896 std::unique_ptr<const TrackParameters> startPar =
unique_clone(cache.m_idmat ? lastidpar.get() : indettrack->perigeeParameters());
898 for (
int i = 0;
i < 2;
i++) {
899 std::unique_ptr<const TrackParameters> tmpelosspar;
901 params1[
Trk::
phi] += muonscatphi;
902 params1[
Trk::
theta] += muonscattheta;
908 std::unique_ptr<const TrackParameters> tmppar1(muonscatpar->associatedSurface().createUniqueTrackParameters(
909 params1[0], params1[1], params1[2], params1[3], params1[4], std::nullopt
921 trajectory.m_fieldprop,
930 calomeots[1].associatedSurface(),
932 trajectory.m_fieldprop
936 if ((tmpelosspar ==
nullptr) || (jac1 == std::nullopt)) {
945 const AmgVector(5) & newpars = tmpelosspar->parameters();
946 std::unique_ptr<const TrackParameters> elosspar2(tmpelosspar->associatedSurface().createUniqueTrackParameters(
947 newpars[0], newpars[1], newpars[2], newpars[3], newqoverpid, std::nullopt
951 elosspar = std::move(tmpelosspar);
954 std::unique_ptr<const TrackParameters> scat2(
m_propagator->propagateParameters(
958 calomeots[0].associatedSurface() :
959 calomeots[2].associatedSurface(),
962 trajectory.m_fieldprop,
972 calomeots[0].associatedSurface() :
973 calomeots[2].associatedSurface(),
977 trajectory.m_fieldprop
981 if ((scat2 ==
nullptr) || (jac2 == std::nullopt)) {
986 for (
int j = 0; j < 5; j++) {
987 for (
int k = 0;
k < 5;
k++) {
989 for (
int l = 0;
l < 5;
l++) {
990 jac3[j][
k] += (*jac2) (j,
l) * (*jac1) (
l,
k);
999 jac4(0, 0) = jac3[0][2];
1000 jac4(1, 1) = jac3[1][3];
1001 jac4(0, 1) = jac3[0][3];
1002 jac4(1, 0) = jac3[1][2];
1004 jac4 = jac4.inverse();
1016 discsurf =
static_cast<const Trk::DiscSurface *
>(&scat2->associatedSurface());
1018 if (cylsurf !=
nullptr) {
1022 if (discsurf !=
nullptr) {
1026 double dphi = jac4(0, 0) * dloc1 + jac4(0, 1) * dloc2;
1027 double dtheta = jac4(1, 0) * dloc1 + jac4(1, 1) * dloc2;
1033 muonscatphi += dphi;
1034 muonscattheta += dtheta;
1037 double idscattheta = idscatpar->parameters()[
Trk::theta] - (scat2->parameters()[
Trk::theta] + dtheta);
1040 firstscatphi = muonscatphi;
1041 secondscatphi = idscatphi;
1042 firstscattheta = muonscattheta;
1043 secondscattheta = idscattheta;
1045 firstscatphi = -idscatphi;
1046 secondscatphi = -muonscatphi;
1047 firstscattheta = -idscattheta;
1048 secondscattheta = -muonscattheta;
1051 if (
i == 1 && cache.m_field_cache.toroidOn() && !firstismuon) {
1053 params2[
Trk::
phi] += idscatphi;
1054 params2[
Trk::
theta] += idscattheta;
1060 firstscatpar=scat2->associatedSurface().createUniqueTrackParameters(
1061 params2[0], params2[1], params2[2], params2[3], params2[4], std::nullopt
1063 idscatpar = firstscatpar.get();
1067 *cache.m_caloEntrance,
1071 if (startPar !=
nullptr) {
1080 rot.col(2) = trackdir;
1083 trans.linear().matrix() << rot;
1084 trans.translation() << startPar->position() - .1 * trackdir;
1085 PlaneSurface curvlinsurf(trans);
1095 if (curvlinpar !=
nullptr) {
1096 startPar.reset(curvlinpar);
1100 firstscatpar = std::move(scat2);
1104 std::unique_ptr<GXFMaterialEffects> firstscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[0]);
1105 std::unique_ptr<GXFMaterialEffects> elossmeff = std::make_unique<GXFMaterialEffects>(calomeots[1]);
1106 std::unique_ptr<GXFMaterialEffects> secondscatmeff = std::make_unique<GXFMaterialEffects>(calomeots[2]);
1108 double pull1 = std::abs(firstscatphi / firstscatmeff->sigmaDeltaPhi());
1109 double pull2 = std::abs(secondscatphi / secondscatmeff->sigmaDeltaPhi());
1112 for (
auto &
i : tmp_matvec) {
1117 firstscatmeff->setScatteringAngles(firstscatphi, firstscattheta);
1118 secondscatmeff->setScatteringAngles(secondscatphi, secondscattheta);
1121 elossmeff->setdelta_p(1000 * (lastscatpar->parameters()[
Trk::qOverP] - newqoverpid));
1123 elossmeff->setdelta_p(1000 * (newqoverpid - firstscatpar->parameters()[
Trk::qOverP]));
1126 elossmeff->setSigmaDeltaE(calomeots[1].energyLoss()->sigmaDeltaE());
1128 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(firstscatmeff), std::move(firstscatpar)), -1);
1129 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(elossmeff), std::move(elosspar)), -1);
1130 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(secondscatmeff), std::move(lastscatpar)), -1);
1133 for (
auto &
i : tmp_matvec) {
1140 if (startPar ==
nullptr) {
1145 (pull1 > 5 || pull2 > 5) &&
1151 bool largegap =
false;
1152 double previousz = 0;
1154 for (itStates2 = beginStates2; itStates2 != endState2; ++itStates2) {
1155 const MaterialEffectsBase *meff = (*itStates2)->materialEffectsOnTrack();
1157 const MeasurementBase *meas = (*itStates2)->measurementOnTrack();
1159 if (meff !=
nullptr) {
1161 const MaterialEffectsOnTrack *mefot{};
1163 mefot =
static_cast<const MaterialEffectsOnTrack *
>(meff);
1165 if ( mefot and mefot->energyLoss() and
1166 std::abs(mefot->energyLoss()->deltaE()) > 250 &&
1167 mefot->energyLoss()->sigmaDeltaE() < 1.e-9
1172 cache.m_extmat =
false;
1174 cache.m_idmat =
false;
1180 !(itStates2 == beginStates2 || itStates2 == beginStates2 + 1) &&
1189 itStates2 == endState2 - 1 &&
1191 tpar->position().perp() > 9000 &&
1192 std::abs(tpar->position().z()) < 13000
1194 std::unique_ptr<const TrackParameters> pseudopar(tpar->clone());
1198 std::unique_ptr<const PseudoMeasurementOnTrack> newpseudo = std::make_unique<const PseudoMeasurementOnTrack>(
1201 pseudopar->associatedSurface()
1204 std::unique_ptr<GXFTrackState> pseudostate = std::make_unique<GXFTrackState>(std::move(newpseudo), std::move(pseudopar));
1211 pseudostate->setMeasurementErrors(
errors);
1212 trajectory.addMeasurementState(std::move(pseudostate));
1218 std::abs(trajectory.trackStates().back()->position().z()) > 20000 &&
1219 std::abs(previousz) < 12000
1225 previousz = trajectory.trackStates().back()->position().z();
1235 (cache.m_field_cache.toroidOn() || cache.m_field_cache.solenoidOn()) ?
muon :
nonInteracting
◆ makePerigee()
Definition at line 4470 of file GlobalChi2Fitter.cxx.
4475 const PerigeeSurface *persurf =
nullptr;
4478 persurf =
static_cast<const PerigeeSurface *
>(¶m.associatedSurface());
4480 if ((persurf !=
nullptr) && (!cache.m_acceleration || persurf->center().perp() > 5)) {
4482 return param.associatedSurface().createUniqueTrackParameters(
4487 if (cache.m_acceleration) {
4491 PerigeeSurface tmppersf;
4492 std::unique_ptr<const TrackParameters> per(
m_extrapolator->extrapolate(
4493 Gaudi::Hive::currentContext(),param, tmppersf,
oppositeMomentum,
false, matEffects));
4495 if (per ==
nullptr) {
4496 ATH_MSG_DEBUG(
"Cannot make Perigee with starting parameters");
4500 if(std::abs(per->position().z())>5000.) {
4501 ATH_MSG_WARNING(
"Pathological perigee well outside of tracking detector!! Returning nullptr");
◆ makeProtoState()
Definition at line 2557 of file GlobalChi2Fitter.cxx.
2569 ) && cache.m_getmaterialfromtrack
2571 if (cache.m_acceleration && trajectory.numberOfHits() == 0) {
2577 const MaterialEffectsOnTrack *meff =
static_cast<const MaterialEffectsOnTrack *
>(tsos->materialEffectsOnTrack());
2579 std::unique_ptr<GXFMaterialEffects> newmeff;
2582 meff->scatteringAngles() or
2583 meff->energyLoss() or
2585 (tsos->trackParameters() ==
nullptr)
2587 newmeff = std::make_unique<GXFMaterialEffects>(*meff);
2591 double sigmascat = std::sqrt(
m_scattool->sigmaSquare(
2593 std::abs(1. / tsos->trackParameters()->parameters()[
Trk::qOverP]),
2606 meff->thicknessInX0(),
2612 newmeff = std::make_unique<GXFMaterialEffects>(newmeot);
2616 (meff->energyLoss() !=
nullptr) &&
2617 meff->energyLoss()->sigmaDeltaE() > 0 &&
2620 ((meff->scatteringAngles() ==
nullptr) || meff->thicknessInX0() == 0)
2623 newmeff->setSigmaDeltaE(meff->energyLoss()->sigmaDeltaE());
2626 (tsos->trackParameters() !=
nullptr) &&
2627 !trajectory.trackStates().empty() &&
2628 ((**trajectory.trackStates().rbegin()).trackParameters() !=
nullptr)
2630 double delta_p = 1000 * (
2631 tsos->trackParameters()->parameters()[
Trk::qOverP] -
2632 (**trajectory.trackStates().rbegin()).trackParameters()->
2636 newmeff->setdelta_p(delta_p);
2640 trajectory.addMaterialState(std::make_unique<GXFTrackState>(std::move(newmeff),
unique_clone(tsos->trackParameters())),
index);
2652 tsos->measurementOnTrack(),
2653 tsos->trackParameters(),
◆ makeProtoStateFromMeasurement()
Definition at line 2660 of file GlobalChi2Fitter.cxx.
2670 if (!measbase->associatedSurface().associatedDetectorElementIdentifier().is_valid()) {
2672 seg =
static_cast<const Segment *
>(measbase);
2679 imax = (
int) seg->numberOfMeasurementBases();
2682 for (
int i = 0;
i <
imax;
i++) {
2683 const MeasurementBase *measbase2 = ((seg !=
nullptr) &&
m_decomposesegments) ? seg->measurement(
i) : measbase;
2685 std::unique_ptr<GXFTrackState> ptsos = std::make_unique<GXFTrackState>(
2686 std::unique_ptr<const MeasurementBase>(measbase2->clone()),
2687 std::unique_ptr<const TrackParameters>(newtrackpar !=
nullptr ? newtrackpar->clone() :
nullptr)
2690 double sinstereo = 0;
2694 Identifier hitid = measbase2->associatedSurface().associatedDetectorElementIdentifier();
2699 hitid = crot->rioOnTrack(0).identify();
2703 bool measphi =
false;
2706 bool rotated =
false;
2722 if (measbase->localParameters().parameterKey() != 1) {
2743 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
2744 errors[0] = std::sqrt(covEigenValueSmall);
2745 sinstereo =
std::sin(covStereoAngle);
2758 const Surface *surf = &measbase2->associatedSurface();
2759 Amg::Vector3D measdir = surf->transform().rotation().col(0);
2761 double dotprod2 = measdir.dot(
Amg::Vector3D(surf->center().x(), surf->center().y(), 0) / surf->center().perp());
2762 if (std::abs(dotprod1) < .5 && std::abs(dotprod2) < .5) {
2770 int param_index = 0;
2777 errors[1] = std::sqrt(
covmat(param_index, param_index));
2782 errors[2] = std::sqrt(
covmat(param_index, param_index));
2787 errors[3] = std::sqrt(
covmat(param_index, param_index));
2792 errors[4] = std::sqrt(
covmat(param_index, param_index));
2797 ATH_MSG_DEBUG(
"PseudoMeasurement, pos=" << measbase2->globalPosition());
2800 ATH_MSG_DEBUG(
"VertexOnTrack, pos=" << measbase2->globalPosition());
2803 ATH_MSG_DEBUG(
"Segment, pos=" << measbase2->globalPosition());
2813 ptsos->setMeasurementErrors(
errors);
2814 ptsos->setSinStereo(sinstereo);
2815 ptsos->setMeasurementType(hittype);
2816 ptsos->setMeasuresPhi(measphi);
2818 if (isoutlier && !cache.m_reintoutl) {
2823 bool ok = trajectory.addMeasurementState(std::move(ptsos),
index);
◆ makeTrack()
Definition at line 7536 of file GlobalChi2Fitter.cxx.
7543 auto trajectory = std::make_unique<Trk::TrackStates>();
7549 GXFTrajectory tmptrajectory(oldtrajectory);
7551 std::unique_ptr<GXFTrackState> perigee_ts =
makeTrackFindPerigee(ctx, cache, oldtrajectory, matEffects);
7553 if (perigee_ts ==
nullptr) {
7557 tmptrajectory.addBasicState(std::move(perigee_ts), cache.m_acceleration ? 0 : tmptrajectory.numberOfUpstreamStates());
7559 trajectory->reserve(tmptrajectory.trackStates().size());
7560 for (
auto & hit : tmptrajectory.trackStates()) {
7565 hit->resetTrackCovariance();
7571 hit->materialEffects()))
7577 auto trackState = hit->trackStateOnSurface();
7578 hit->resetTrackCovariance();
7579 trajectory->emplace_back(trackState.release());
7582 auto qual = std::make_unique<FitQuality>(tmptrajectory.chi2(), tmptrajectory.nDOF());
7593 if (matEffects ==
electron && tmptrajectory.hasKink()) {
7598 if (tmptrajectory.m_straightline) {
7602 std::unique_ptr<Track>
rv = std::make_unique<Track>(
info, std::move(trajectory), std::move(
qual));
7612 std::unique_ptr<TrackSummary>
ts = std::make_unique<TrackSummary>();
7621 std::optional<TrackHoleCount> hole_count;
7646 if (hole_count.has_value()) {
7659 rv->setTrackSummary(std::move(
ts));
◆ makeTrackFillDerivativeMatrix()
void Trk::GlobalChi2Fitter::makeTrackFillDerivativeMatrix |
( |
Cache & |
cache, |
|
|
GXFTrajectory & |
oldtrajectory |
|
) |
| |
|
staticprivate |
Definition at line 6996 of file GlobalChi2Fitter.cxx.
7000 Amg::MatrixX & derivs = oldtrajectory.weightedResidualDerivatives();
7004 for (
auto & hit : oldtrajectory.trackStates()) {
7005 if (
const auto *pMeas{hit->measurement()};
7011 nrealmeas += hit->numberOfMeasuredParameters();
7014 cache.m_derivmat.resize(nrealmeas, oldtrajectory.numberOfFitParameters());
7015 cache.m_derivmat.setZero();
7018 int nperpars = oldtrajectory.numberOfPerigeeParameters();
7019 int nscat = oldtrajectory.numberOfScatterers();
7020 for (
auto & hit : oldtrajectory.trackStates()) {
7021 if (
const auto *pMeas{hit->measurement()};
7027 for (
int i = measindex;
i < measindex + hit->numberOfMeasuredParameters();
i++) {
7028 for (
int j = 0; j < oldtrajectory.numberOfFitParameters(); j++) {
7029 cache.m_derivmat(
i, j) = derivs(measindex2, j) *
errors[measindex2];
7030 if ((j == 4 && !oldtrajectory.m_straightline) || j >= nperpars + 2 * nscat) {
7031 cache.m_derivmat(
i, j) *= 1000;
7038 measindex += hit->numberOfMeasuredParameters();
7039 }
else if (hit->materialEffects() ==
nullptr) {
7040 measindex2 += hit->numberOfMeasuredParameters();
◆ makeTrackFindPerigee()
◆ makeTrackFindPerigeeParameters()
Definition at line 7045 of file GlobalChi2Fitter.cxx.
7051 GXFTrackState *firstmeasstate =
nullptr, *lastmeasstate =
nullptr;
7052 std::tie(firstmeasstate, lastmeasstate) = oldtrajectory.findFirstLastMeasurement();
7053 std::unique_ptr<const TrackParameters> per(
nullptr);
7056 std::unique_ptr<const TrackParameters> prevpar(
7057 firstmeasstate->trackParameters() !=
nullptr ?
7058 firstmeasstate->trackParameters()->clone() :
7061 std::vector<std::pair<const Layer *, const Layer *>> & upstreamlayers = oldtrajectory.upstreamMaterialLayers();
7065 if (prevpar ==
nullptr) {
7070 const Layer *
layer = layer1 !=
nullptr ? layer1 : layer2;
7072 DistanceSolution distsol =
layer->surfaceRepresentation().straightLineDistanceEstimate(
7073 prevpar->position(), prevpar->momentum().unit()
7077 if (distsol.numberOfSolutions() == 2) {
7082 if (distsol.first() * distsol.second() < 0 && !
first) {
7091 std::unique_ptr<const TrackParameters> layerpar(
7095 layer->surfaceRepresentation(),
7098 oldtrajectory.m_fieldprop,
7103 if (layerpar ==
nullptr) {
7107 if (
layer->surfaceRepresentation().bounds().inside(layerpar->localPosition())) {
7111 prevpar = std::move(layerpar);
7115 const Layer *startlayer = firstmeasstate->trackParameters()->associatedSurface().associatedLayer();
7117 if ((startlayer !=
nullptr) && (startlayer->layerMaterialProperties() !=
nullptr)) {
7118 double startfactor = startlayer->layerMaterialProperties()->alongPostFactor();
7119 const Surface & discsurf = startlayer->surfaceRepresentation();
7122 startfactor = startlayer->layerMaterialProperties()->oppositePostFactor();
7124 if (startfactor > 0.5) {
7125 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7126 firstmeasstate->trackParameters(), *startlayer,
oppositeMomentum, matEffects
7129 if (updatedpar !=
nullptr) {
7130 firstmeasstate->setTrackParameters(std::move(updatedpar));
7139 const Layer *endlayer = lastmeasstate->trackParameters()->associatedSurface().associatedLayer();
7141 if ((endlayer !=
nullptr) && (endlayer->layerMaterialProperties() !=
nullptr)) {
7142 double endfactor = endlayer->layerMaterialProperties()->alongPreFactor();
7143 const Surface & discsurf = endlayer->surfaceRepresentation();
7146 endfactor = endlayer->layerMaterialProperties()->oppositePreFactor();
7149 if (endfactor > 0.5) {
7150 std::unique_ptr<const TrackParameters> updatedpar =
m_matupdator->update(
7151 lastmeasstate->trackParameters(), *endlayer,
alongMomentum, matEffects
7154 if (updatedpar !=
nullptr) {
7155 lastmeasstate->setTrackParameters(std::move(updatedpar));
7160 if (prevpar !=
nullptr) {
7167 oldtrajectory.m_fieldprop,
7172 if (per ==
nullptr) {
7173 ATH_MSG_DEBUG(
"Failed to extrapolate to perigee, returning 0");
7178 }
else if (cache.m_acceleration && (firstmeasstate->trackParameters() !=
nullptr)) {
7180 *firstmeasstate->trackParameters(),
7186 per.reset(oldtrajectory.referenceParameters()->clone());
◆ myfit()
Definition at line 4508 of file GlobalChi2Fitter.cxx.
4516 ATH_MSG_DEBUG(
"--> entering GlobalChi2Fitter::myfit_helper");
4520 if (!trajectory.m_straightline) {
4521 if (trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
4522 trajectory.m_straightline = !cache.m_field_cache.solenoidOn();
4523 }
else if ((trajectory.prefit() == 0) && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == 0) {
4524 trajectory.m_straightline = !cache.m_field_cache.toroidOn();
4526 trajectory.m_straightline = (!cache.m_field_cache.solenoidOn() && !cache.m_field_cache.toroidOn());
4531 cache.m_lastiter = 0;
4535 if (trajectory.numberOfPerigeeParameters() == -1) {
4536 cache.incrementFitStatus(
S_FITS);
4537 if (trajectory.m_straightline) {
4538 trajectory.setNumberOfPerigeeParameters(4);
4540 trajectory.setNumberOfPerigeeParameters(5);
4544 if (trajectory.nDOF() < 0) {
4549 cache.m_phiweight.clear();
4550 cache.m_firstmeasurement.clear();
4551 cache.m_lastmeasurement.clear();
4554 ATH_MSG_WARNING(
"Attempt to apply material corrections with q/p=0, reject track");
4558 if (matEffects ==
Trk::electron && trajectory.m_straightline) {
4564 trajectory.setMass(
mass);
4566 ATH_MSG_DEBUG(
"start param: " << param <<
" pos: " << param.position() <<
" pt: " << param.pT());
4568 std::unique_ptr<const TrackParameters> per =
makePerigee(cache, param, matEffects);
4570 if (!cache.m_acceleration && (per ==
nullptr)) {
4580 cache.m_acceleration &&
4581 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits() &&
4584 ATH_MSG_WARNING(
"Tracking Geometry Service and/or Material Updator Tool not configured");
4587 cache.m_fastmat =
false;
4592 !cache.m_acceleration ||
4593 trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() != trajectory.numberOfHits()
4595 addMaterial(ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4598 ctx, cache, trajectory, per !=
nullptr ? per.get() : ¶m, matEffects);
4602 if (cache.m_acceleration && (trajectory.referenceParameters() ==
nullptr) && (per ==
nullptr)) {
4605 if (trajectory.numberOfScatterers() >= 2) {
4606 GXFTrackState *scatstate =
nullptr;
4607 GXFTrackState *scatstate2 =
nullptr;
4610 for (
const auto & state : trajectory.trackStates()) {
4613 scatindex == trajectory.numberOfScatterers() / 2 ||
4614 state->materialEffects()->deltaE() == 0
4616 scatstate2 = state.get();
4621 scatstate = state.get();
4628 if ((scatstate ==
nullptr) || (scatstate2 ==
nullptr)) {
4629 throw std::logic_error(
"Invalid scatterer");
4632 vertex = .49 * (scatstate->position() + scatstate2->position());
4634 int nstates = (
int) trajectory.trackStates().size();
4636 trajectory.trackStates()[nstates / 2 - 1]->position() +
4637 trajectory.trackStates()[nstates / 2]->position()
4641 PerigeeSurface persurf(
vertex);
4642 std::unique_ptr<const TrackParameters> nearestpar;
4643 double mindist = 99999;
4644 std::vector < GXFTrackState * >mymatvec;
4646 for (
auto &
it : trajectory.trackStates()) {
4647 if ((*it).trackParameters() ==
nullptr) {
4652 .straightLineDistanceEstimate(
4653 (*it).trackParameters()->position(),
4654 (*it).trackParameters()->momentum().unit())
4658 (cache.m_caloEntrance ==
nullptr) ||
4659 cache.m_caloEntrance->inside((*it).trackParameters()->position())
4663 (((*it).measurement() !=
nullptr) && insideid) || (
4664 ((*it).materialEffects() !=
nullptr) &&
4666 (*it).materialEffects()->deltaE() == 0 ||
4667 ((*it).materialEffects()->sigmaDeltaPhi() == 0 &&
4669 (*it).materialEffects()->deltaPhi() != 0
4673 double dist = ((*it).trackParameters()->position() -
vertex).
perp();
4674 if (dist < mindist) {
4682 if (((*it).materialEffects() !=
nullptr) &&
distance > 0) {
4683 mymatvec.push_back(
it.get());
4687 if (nearestpar ==
nullptr) {
4691 for (
auto & state : mymatvec) {
4693 const Surface &matsurf = state->associatedSurface();
4694 DistanceSolution distsol = matsurf.straightLineDistanceEstimate(
4695 nearestpar->position(), nearestpar->momentum().unit());
4699 if (
distance < 0 && distsol.numberOfSolutions() > 0) {
4703 std::unique_ptr<const TrackParameters> tmppar(
m_propagator->propagateParameters(
4709 trajectory.m_fieldprop,
4713 if (tmppar ==
nullptr) {
4721 trajectory.m_fieldprop,
4725 if (tmppar ==
nullptr) {
4737 if (state->materialEffects()->sigmaDeltaE() > 0) {
4738 newpars[
Trk::qOverP] += .001 * state->materialEffects()->delta_p();
4741 double de = std::abs(state->materialEffects()->deltaE());
4743 double newp2 = oldp * oldp - 2 * de * std::sqrt(
mass *
mass + oldp * oldp) + de * de;
4749 nearestpar = tmppar->associatedSurface().createUniqueTrackParameters(
4750 newpars[0], newpars[1], newpars[2], newpars[3], newpars[4], std::nullopt
4754 std::unique_ptr<Trk::TrackParameters> tmpPars(
m_propagator->propagateParameters(
4760 trajectory.m_fieldprop,
4765 if (tmpPars !=
nullptr) {
4766 per.reset(
static_cast < const Perigee *
>(tmpPars.release()));
4771 double oldp = 1. / std::abs(per->parameters()[
Trk::qOverP]);
4772 double toteloss = std::abs(trajectory.totalEnergyLoss());
4773 double newp = std::sqrt(oldp * oldp + 2 * toteloss * std::sqrt(oldp * oldp +
mass *
mass) + toteloss * toteloss);
4777 per = per->associatedSurface().createUniqueTrackParameters(
4782 if (per ==
nullptr) {
4790 PerigeeSurface persurf2(per->position());
4791 per = persurf2.createUniqueTrackParameters(
4799 }
else if (per ==
nullptr) {
4803 if ((per ==
nullptr) && (trajectory.referenceParameters() ==
nullptr)) {
4811 if (trajectory.m_straightline && (per !=
nullptr)) {
4812 if (trajectory.numberOfPerigeeParameters() == -1) {
4813 trajectory.setNumberOfPerigeeParameters(4);
4817 per = per->associatedSurface().createUniqueTrackParameters(
4820 }
else if (trajectory.numberOfPerigeeParameters() == -1) {
4821 trajectory.setNumberOfPerigeeParameters(5);
4824 if (per !=
nullptr) {
4825 trajectory.setReferenceParameters(std::move(per));
4828 int nfitpar = trajectory.numberOfFitParameters();
4829 int nperpars = trajectory.numberOfPerigeeParameters();
4830 int nscat = trajectory.numberOfScatterers();
4831 int nbrem = trajectory.numberOfBrems();
4834 Eigen::MatrixXd a_inv;
4835 a.resize(nfitpar, nfitpar);
4840 derivPool.setZero();
4842 for (std::unique_ptr<GXFTrackState> & state : trajectory.trackStates()) {
4843 if (state->materialEffects() !=
nullptr) {
4846 state->setDerivatives(derivPool);
4849 bool doderiv =
true;
4850 int tmpminiter = cache.m_miniter;
4853 cache.m_lastiter =
it;
4859 cache.m_miniter = tmpminiter;
4863 if (!trajectory.converged()) {
4864 cache.m_fittercode =
4874 cache.m_miniter = tmpminiter;
4878 int nhits = trajectory.numberOfHits();
4879 int ntrthits = trajectory.numberOfTRTHits();
4880 int nsihits = trajectory.numberOfSiliconHits();
4881 double redchi2 = (trajectory.nDOF() > 0) ? trajectory.chi2() / trajectory.nDOF() : 0;
4882 double prevredchi2 = (trajectory.nDOF() > 0) ? trajectory.prevchi2() / trajectory.nDOF() : 0;
4891 (redchi2 < prevredchi2 &&
4892 (redchi2 > prevredchi2 - 1 || redchi2 < 2)) ||
4893 nsihits + ntrthits == nhits
4898 if (
it != 1 || nsihits != 0 || trajectory.nDOF() <= 0 || trajectory.chi2() / trajectory.nDOF() <= 3) {
4903 cache.m_miniter = tmpminiter;
4910 int ntrtprechits = trajectory.numberOfTRTPrecHits();
4911 int ntrttubehits = trajectory.numberOfTRTTubeHits();
4913 if (ntrtprechits+ntrttubehits) {
4914 phf =
float(ntrtprechits)/
float(ntrtprechits+ntrttubehits);
4916 if (phf<m_minphfcut && it>=3) {
4917 if ((ntrtprechits+ntrttubehits)>=15) {
4922 <<
" | nTRTPrecHits = " << ntrtprechits
4923 <<
" | nTRTTubeHits = " << ntrttubehits
4924 <<
" | nOutliers = "
4925 << trajectory.numberOfOutliers());
4927 if (!trajectory.converged()) {
4933 cache.m_miniter = tmpminiter;
4942 cache.m_miniter = tmpminiter;
4944 if (trajectory.prefit() == 0) {
4948 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
4949 if (lltOfW.info() == Eigen::Success) {
4953 int ncols =
a.cols();
4954 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
4955 a_inv = lltOfW.solve(weightInvAMG);
4964 GXFTrajectory *finaltrajectory = &trajectory;
4966 (runOutlier || cache.m_sirecal) &&
4967 trajectory.numberOfSiliconHits() == trajectory.numberOfHits()
4974 if (traj != &trajectory) {
4979 finaltrajectory = traj;
4986 if (!cache.m_acceleration && (finaltrajectory->prefit() == 0)) {
4987 if (nperpars == 5) {
4988 for (
int i = 0;
i <
a.cols();
i++) {
4989 a_inv(4,
i) *= .001;
4990 a_inv(
i, 4) *= .001;
4994 int scatterPos = nperpars + 2 * nscat;
4995 for (
int bremno = 0; bremno < nbrem; bremno++, scatterPos++) {
4996 for (
int i = 0;
i <
a.cols();
i++) {
4997 a_inv(scatterPos,
i) *= .001;
4998 a_inv(
i, scatterPos) *= .001;
5004 int nperparams = finaltrajectory->numberOfPerigeeParameters();
5005 for (
int i = 0;
i < nperparams;
i++) {
5006 for (
int j = 0; j < nperparams; j++) {
5007 (errmat) (j,
i) = a_inv(j,
i);
5011 if (trajectory.m_straightline) {
5012 (errmat) (4, 4) = 1
e-20;
5015 const AmgVector(5) & perpars = finaltrajectory->referenceParameters()->parameters();
5016 std::unique_ptr<const TrackParameters> measper(
5017 finaltrajectory->referenceParameters()->associatedSurface().createUniqueTrackParameters(
5018 perpars[0], perpars[1], perpars[2], perpars[3], perpars[4], std::move(errmat)
5022 finaltrajectory->setReferenceParameters(std::move(measper));
5024 cache.m_fullcovmat = a_inv;
5028 std::unique_ptr<Track>
track =
nullptr;
5030 if (finaltrajectory->prefit() > 0) {
5031 if (finaltrajectory != &trajectory) {
5033 delete finaltrajectory;
5038 if (finaltrajectory->numberOfOutliers() <=
m_maxoutliers || !runOutlier) {
5045 double cut = (finaltrajectory->numberOfSiliconHits() ==
5046 finaltrajectory->numberOfHits())
5052 (track !=
nullptr) && (
5053 track->fitQuality()->numberDoF() != 0 &&
5054 track->fitQuality()->chiSquared() /
track->fitQuality()->numberDoF() >
cut
5057 track.reset(
nullptr);
5061 if (track ==
nullptr) {
5065 if (finaltrajectory != &trajectory) {
5066 delete finaltrajectory;
5069 return track.release();
◆ myfit_helper()
◆ numericalDerivatives()
Definition at line 8301 of file GlobalChi2Fitter.cxx.
8308 ParamDefsAccessor paraccessor;
8316 std::optional<TransportJacobian> jac = std::make_optional<TransportJacobian>(J);
8319 0.01, 0.01, 0.00001, 0.00001, 0.000000001
8326 const Surface & previousSurface = tmpprevpar->associatedSurface();
8330 for (
int i = 0;
i < 5;
i++) {
8333 if (thisdiscsurf &&
i == 1) {
8337 vecpluseps[paraccessor.pardef[
i]] += eps[
i];
8338 vecminuseps[paraccessor.pardef[
i]] -= eps[
i];
8339 if (
i == 0 && thiscylsurf) {
8341 }
else if (
i == 1 && thisdiscsurf) {
8347 std::unique_ptr<const TrackParameters> parpluseps(
8348 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8357 std::unique_ptr<const TrackParameters> parminuseps(
8358 tmpprevpar->associatedSurface().createUniqueTrackParameters(
8368 std::unique_ptr<const TrackParameters> newparpluseps(
8379 std::unique_ptr<const TrackParameters> newparminuseps(
8394 if (newparpluseps ==
nullptr) {
8406 if (newparminuseps ==
nullptr) {
8418 if ((newparpluseps ==
nullptr) || (newparminuseps ==
nullptr)) {
8422 for (
int j = 0; j < 5; j++) {
8423 double diff = newparpluseps->parameters()[paraccessor.pardef[j]] -
8424 newparminuseps->parameters()[paraccessor.pardef[j]];
8426 if (j == 0 && cylsurf) {
8428 }
else if (j == 1 && discsurf) {
8432 (*jac) (j,
i) =
diff / (2 * eps[
i]);
◆ processTrkVolume()
Definition at line 2833 of file GlobalChi2Fitter.cxx.
2837 if (tvol ==
nullptr) {
2844 if (confinedLayers !=
nullptr) {
2848 if (
layer !=
nullptr) {
2853 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2857 const CylinderLayer *cyllay =
nullptr;
2859 cyllay =
static_cast<const CylinderLayer *
>(
layer);
2862 const DiscLayer *disclay =
nullptr;
2864 disclay =
static_cast<const DiscLayer *
>(
layer);
2867 if (disclay !=
nullptr) {
2868 if (disclay->center().z() < 0) {
2869 cache.m_negdiscs.push_back(disclay);
2871 cache.m_posdiscs.push_back(disclay);
2873 }
else if (cyllay !=
nullptr) {
2874 cache.m_barrelcylinders.push_back(cyllay);
2884 for (
size_t ib = 0 ;
ib < bsurf.size(); ++
ib) {
2885 const Layer *
layer = bsurf[
ib]->surfaceRepresentation().materialLayer();
2887 if (
layer ==
nullptr)
continue;
2891 if ((layIndex.
value() == 0) || (
layer->layerMaterialProperties() ==
nullptr)) {
2895 const CylinderSurface *cylsurf =
nullptr;
2898 cylsurf =
static_cast<const CylinderSurface *
>(&
layer->surfaceRepresentation());
2900 const DiscSurface *discsurf =
nullptr;
2903 discsurf =
static_cast<const DiscSurface *
>(&
layer->surfaceRepresentation());
2905 if (discsurf !=
nullptr) {
2907 discsurf->center().z() < 0 &&
2908 std::find(cache.m_negdiscs.begin(), cache.m_negdiscs.end(),
layer) == cache.m_negdiscs.end()
2910 cache.m_negdiscs.push_back(
layer);
2912 discsurf->center().z() > 0 &&
2913 std::find(cache.m_posdiscs.begin(), cache.m_posdiscs.end(),
layer) == cache.m_posdiscs.end()
2915 cache.m_posdiscs.push_back(
layer);
2918 (cylsurf !=
nullptr) &&
2919 std::find(cache.m_barrelcylinders.begin(), cache.m_barrelcylinders.end(),
layer) == cache.m_barrelcylinders.end()
2921 cache.m_barrelcylinders.push_back(
layer);
2924 if ((cylsurf ==
nullptr) && (discsurf ==
nullptr)) {
2931 if (confinedVolumes !=
nullptr) {
2932 for (
const auto & volume : confinedVolumes->arrayObjects()) {
2933 if (volume !=
nullptr) {
◆ retrieveTrackingGeometry()
◆ runIteration()
Definition at line 6000 of file GlobalChi2Fitter.cxx.
6010 const int nDOFold = trajectory.nDOF();
6011 const double oldChi2 = trajectory.chi2();
6012 const double oldRedChi2 = nDOFold > 0 ? oldChi2 / nDOFold : 0;
6014 if (cache.m_phiweight.empty()) {
6015 cache.m_phiweight.assign(trajectory.trackStates().size(), 1);
6034 int bremno_maxbrempull = 0;
6035 GXFTrackState* state_maxbrempull =
nullptr;
6050 if ((state_maxbrempull !=
nullptr) && trajectory.converged()) {
6051 trajectory.setConverged(
false);
6052 trajectory.setChi2(1e15);
6059 const int nDOFnew = trajectory.nDOF();
6060 const double newChi2 = trajectory.chi2();
6061 const double newRedChi2 = nDOFnew > 0 ? newChi2 / nDOFnew : 0;
6063 ATH_MSG_DEBUG(
"old chi2: " << oldChi2 <<
"/" << nDOFold <<
"=" << oldRedChi2 <<
6064 ", new chi2: " << newChi2 <<
"/" << nDOFnew <<
"=" << newRedChi2);
6066 if (trajectory.prefit() > 0 && trajectory.converged()) {
6075 if (cache.m_firstmeasurement.empty()) {
6079 if (
a.cols() != trajectory.numberOfFitParameters()) {
6099 if (doDeriv || weightChanged) {
6110 if (trajectory.prefit() == 0) {
6111 if (trajectory.converged()) {
6112 const int nSiHits = trajectory.numberOfSiliconHits();
6113 const int nTrtHits = trajectory.numberOfTRTHits();
6114 const int nHits = trajectory.numberOfHits();
6116 if (
nSiHits + nTrtHits != nHits) {
6123 (newRedChi2 < 2 || (newRedChi2 < oldRedChi2 && newRedChi2 > oldRedChi2 - .5))
◆ runTrackCleanerSilicon()
- Warning
- This method has some unclear memory ownership mechanics that might not correspond fully with the model described at the beginning of the file. Be aware!
Definition at line 6517 of file GlobalChi2Fitter.cxx.
6526 bool trackok =
false;
6527 GXFTrajectory *oldtrajectory = &trajectory;
6528 std::unique_ptr < GXFTrajectory > cleanup_oldtrajectory;
6529 GXFTrajectory *newtrajectory =
nullptr;
6530 std::unique_ptr < GXFTrajectory > cleanup_newtrajectory;
6536 while (!trackok && oldtrajectory->nDOF() > 0) {
6538 std::vector<std::unique_ptr<GXFTrackState>> &
states = oldtrajectory->trackStates();
6541 Amg::MatrixX & weightderiv = oldtrajectory->weightedResidualDerivatives();
6542 int nfitpars = oldtrajectory->numberOfFitParameters();
6543 int nhits = oldtrajectory->numberOfHits();
6544 int nsihits = oldtrajectory->numberOfSiliconHits();
6546 if (nhits != nsihits) {
6550 double maxsipull = -1;
6552 int hitno_maxsipull = -1;
6553 int measno_maxsipull = -1;
6554 int stateno_maxsipull = 0;
6555 GXFTrackState *state_maxsipull =
nullptr;
6560 int noutl = oldtrajectory->numberOfOutliers();
6566 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6567 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6573 double *
errors = state->measurementErrors();
6575 const Amg::MatrixX & hitcov = state->measurement()->localCovariance();
6576 double sinstereo = state->sinStereo();
6577 double cosstereo = (sinstereo == 0) ? 1 : std::sqrt(1 - sinstereo * sinstereo);
6578 double weight1 = -1;
6580 if (hitcov(0, 0) > trackcov(0, 0)) {
6581 if (sinstereo == 0) {
6585 trackcov(0, 0) * cosstereo * cosstereo + 2 *
6586 trackcov(1, 0) * cosstereo * sinstereo + trackcov(1, 1) * sinstereo * sinstereo
6597 double sipull1 = weight1 > 0 ? std::abs(
res[measno] / std::sqrt(weight1)) : -1;
6600 std::abs(
res[measno + 1] / std::sqrt(weight2)) :
6603 sipull1 =
std::max(sipull1, sipull2);
6605 if (sipull1 > maxsipull) {
6606 maxsipull = sipull1;
6607 measno_maxsipull = measno;
6608 state_maxsipull = state.get();
6609 stateno_maxsipull = stateno;
6610 hitno_maxsipull = hitno;
6621 measno += state->numberOfMeasuredParameters();
6625 double maxpull = maxsipull;
6627 ATH_MSG_DEBUG(
" maxsipull: " << maxsipull <<
" hitno_maxsipull: " <<
6628 hitno_maxsipull <<
" n3sigma: " << n3sigma <<
" cut: " <<
cut <<
" cut2: " << cut2);
6637 oldtrajectory->chi2() / oldtrajectory->nDOF() > .25 *
m_chi2cut
6639 state_maxsipull = oldtrajectory->trackStates()[stateno_maxsipull].get();
6640 const PrepRawData *prd{};
6642 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6643 prd = rot->prepRawData();
6645 std::unique_ptr < const RIO_OnTrack > broadrot;
6646 double *olderror = state_maxsipull->measurementErrors();
6648 const TrackParameters *trackpar_maxsipull = state_maxsipull->trackParameters();
6650 Amg::VectorX parameterVector = trackpar_maxsipull->parameters();
6651 std::unique_ptr<const TrackParameters> trackparForCorrect(
6652 trackpar_maxsipull->associatedSurface().createUniqueTrackParameters(
6658 state_maxsipull->hasTrackCovariance()
6660 state_maxsipull->trackCovariance())
6664 newerror[0] = newerror[1] = newerror[2] = newerror[3] = newerror[4] = -1;
6665 double newpull = -1;
6666 double newpull1 = -1;
6667 double newpull2 = -1;
6668 double newres1 = -1;
6669 double newres2 = -1;
6670 double newsinstereo = 0;
6674 !state_maxsipull->isRecalibrated() &&
6676 oldtrajectory->chi2() / trajectory.nDOF() > .3 *
m_chi2cut &&
6685 if (state_maxsipull->sinStereo() != 0) {
6686 const auto [covEigenValueSmall, covStereoAngle] = principalComponentAnalysis2x2(
covmat);
6687 newerror[0] = std::sqrt(covEigenValueSmall);
6688 newsinstereo =
std::sin(covStereoAngle);
6690 newerror[0] = std::sqrt(
covmat(0, 0));
6693 double cosstereo = (newsinstereo == 0) ? 1. : std::sqrt(1 - newsinstereo * newsinstereo);
6695 if (cosstereo != 1.) {
6697 cosstereo * (broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX]) +
6698 newsinstereo * (broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY])
6701 newres1 = broadrot->localParameters()[
Trk::locX] - trackpar_maxsipull->parameters()[
Trk::locX];
6704 if (newerror[0] == 0.0) {
6705 ATH_MSG_WARNING(
"Measurement error is zero or negative, treating as outlier");
6708 newpull1 = std::abs(newres1 / newerror[0]);
6712 newerror[1] = std::sqrt(
covmat(1, 1));
6713 newres2 = broadrot->localParameters()[
Trk::locY] - trackpar_maxsipull->parameters()[
Trk::locY];
6714 newpull2 = std::abs(newres2 / newerror[1]);
6717 newpull =
std::max(newpull1, newpull2);
6723 (newerror[0] > 1.5 * olderror[0] || newerror[1] > 1.5 * std::abs(olderror[1]))
6725 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6726 throw std::runtime_error(
6727 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6732 newtrajectory = oldtrajectory;
6734 if (
a.cols() != nfitpars) {
6738 double oldres1 =
res[measno_maxsipull];
6739 res[measno_maxsipull] = newres1;
6740 err[measno_maxsipull] = newerror[0];
6742 for (
int i = 0;
i < nfitpars;
i++) {
6743 if (weightderiv(measno_maxsipull,
i) == 0) {
6747 b[
i] -= weightderiv(measno_maxsipull,
i) * (oldres1 / olderror[0] - (newres1 * olderror[0]) / (newerror[0] * newerror[0]));
6749 for (
int j =
i; j < nfitpars; j++) {
6753 weightderiv(measno_maxsipull,
i) *
6754 weightderiv(measno_maxsipull, j) *
6755 ((olderror[0] * olderror[0]) / (newerror[0] * newerror[0]) - 1)
6759 weightderiv(measno_maxsipull,
i) *= olderror[0] / newerror[0];
6763 double oldres2 =
res[measno_maxsipull + 1];
6764 res[measno_maxsipull + 1] = newres2;
6765 err[measno_maxsipull + 1] = newerror[1];
6767 for (
int i = 0;
i < nfitpars;
i++) {
6768 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6772 b[
i] -= weightderiv(measno_maxsipull + 1,
i) * (oldres2 / olderror[1] - (newres2 * olderror[1]) / (newerror[1] * newerror[1]));
6774 for (
int j =
i; j < nfitpars; j++) {
6778 weightderiv(measno_maxsipull + 1,
i) *
6779 weightderiv(measno_maxsipull + 1, j) *
6780 ((olderror[1] * olderror[1]) / (newerror[1] * newerror[1]) - 1)
6785 weightderiv(measno_maxsipull + 1,
i) *= olderror[1] / newerror[1];
6790 "Recovering outlier, hitno=" << hitno_maxsipull <<
" measno=" <<
6791 measno_maxsipull <<
" pull=" << maxsipull <<
" olderror_0=" <<
6792 olderror[0] <<
" newerror_0=" << newerror[0] <<
" olderror_1=" <<
6793 olderror[1] <<
" newerror_1=" << newerror[1]
6796 state_maxsipull->setMeasurement(std::move(broadrot));
6797 state_maxsipull->setSinStereo(newsinstereo);
6798 state_maxsipull->setMeasurementErrors(newerror);
6802 ((n3sigma < 2 && maxsipull > cut2 && maxsipull <
cut) || n3sigma > 1) &&
6803 (oldtrajectory->chi2() / oldtrajectory->nDOF() > .3 *
m_chi2cut || noutl > 1)
6807 (oldtrajectory->nDOF() > 1 || hittype_maxsipull ==
TrackState::SCT) &&
6812 "Removing outlier, hitno=" << hitno_maxsipull <<
", measno=" <<
6813 measno_maxsipull <<
" pull=" << maxsipull
6820 cleanup_newtrajectory = std::make_unique<GXFTrajectory>(*oldtrajectory);
6821 newtrajectory = cleanup_newtrajectory.get();
6823 if (newa.cols() != nfitpars) {
6828 Amg::MatrixX & newweightderiv = newtrajectory->weightedResidualDerivatives();
6829 if ((measno_maxsipull < 0) or(measno_maxsipull >= (
int)
res.size())) {
6830 throw std::runtime_error(
6831 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6835 double oldres1 =
res[measno_maxsipull];
6836 newres[measno_maxsipull] = 0;
6838 for (
int i = 0;
i < nfitpars;
i++) {
6839 if (weightderiv(measno_maxsipull,
i) == 0) {
6843 newb[
i] -= weightderiv(measno_maxsipull,
i) * oldres1 / olderror[0];
6845 for (
int j =
i; j < nfitpars; j++) {
6849 weightderiv(measno_maxsipull,
i) *
6850 weightderiv(measno_maxsipull, j)
6854 newweightderiv(measno_maxsipull,
i) = 0;
6858 double oldres2 =
res[measno_maxsipull + 1];
6859 newres[measno_maxsipull + 1] = 0;
6861 for (
int i = 0;
i < nfitpars;
i++) {
6862 if (weightderiv(measno_maxsipull + 1,
i) == 0) {
6866 newb[
i] -= weightderiv(measno_maxsipull + 1,
i) * oldres2 / olderror[1];
6868 for (
int j =
i; j < nfitpars; j++) {
6869 if (weightderiv(measno_maxsipull + 1, j) == 0) {
6876 weightderiv(measno_maxsipull + 1,
i) *
6877 weightderiv(measno_maxsipull + 1, j)
6881 newweightderiv(measno_maxsipull + 1,
i) = 0;
6885 newtrajectory->setOutlier(stateno_maxsipull);
6891 newtrajectory->setConverged(
false);
6907 if (!newtrajectory->converged()) {
6909 ctx, cache, *newtrajectory,
it, *newap, *newbp, lu_m, doderiv);
6916 if (!newtrajectory->converged()) {
6925 double oldchi2 = oldtrajectory->chi2() / oldtrajectory->nDOF();
6926 double newchi2 = (newtrajectory->nDOF() > 0) ? newtrajectory->chi2() / newtrajectory->nDOF() : 0;
6929 if (newtrajectory->nDOF() != oldtrajectory->nDOF() && maxsipull > cut2) {
6930 mindiff = (oldchi2 > .33 *
m_chi2cut || noutl > 0) ? .8 : 1.;
6932 if (noutl == 0 && maxsipull <
cut - .5 && oldchi2 < .5 *
m_chi2cut) {
6937 if (newchi2 > oldchi2 || (newchi2 > oldchi2 - mindiff && newchi2 > .33 * oldchi2)) {
6938 ATH_MSG_DEBUG(
"Outlier not confirmed, keeping old trajectory");
6946 (void)cleanup_oldtrajectory.release();
6947 return oldtrajectory;
6949 if (oldtrajectory != newtrajectory) {
6950 cleanup_oldtrajectory = std::move(cleanup_newtrajectory);
6951 oldtrajectory = newtrajectory;
6958 Eigen::LLT < Eigen::MatrixXd > lltOfW(
a);
6959 if (lltOfW.info() == Eigen::Success) {
6963 int ncols =
a.cols();
6964 Amg::MatrixX weightInvAMG = Amg::MatrixX::Identity(ncols, ncols);
6965 fullcov = lltOfW.solve(weightInvAMG);
6983 oldtrajectory->nDOF() > 0 &&
6984 oldtrajectory->chi2() / oldtrajectory->nDOF() >
m_chi2cut &&
6992 (void)cleanup_oldtrajectory.release();
6993 return oldtrajectory;
◆ runTrackCleanerTRT()
Definition at line 6347 of file GlobalChi2Fitter.cxx.
6360 if (
it == 1 && trajectory.numberOfSiliconHits() + trajectory.numberOfTRTHits() == trajectory.numberOfHits()) {
6364 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6367 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6368 int nfitpars = trajectory.numberOfFitParameters();
6370 if (
a.cols() != nfitpars) {
6374 int nperpars = trajectory.numberOfPerigeeParameters();
6375 int nscats = trajectory.numberOfScatterers();
6378 bool outlierremoved =
false;
6379 bool hitrecalibrated =
false;
6381 for (
int stateno = 0; stateno < (
int)
states.size(); stateno++) {
6382 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6390 std::abs(state->trackParameters()->parameters()[
Trk::driftRadius]) > 1.05 * state->associatedSurface().bounds().r()
6394 trajectory.setOutlier(stateno);
6395 outlierremoved =
true;
6397 double *
errors = state->measurementErrors();
6398 double olderror =
errors[0];
6400 trajectory.updateTRTHitCount(stateno, olderror);
6402 for (
int i = 0;
i < nfitpars;
i++) {
6403 if (weightderiv(measno,
i) == 0) {
6407 b[
i] -=
res[measno] * weightderiv(measno,
i) / olderror;
6409 for (
int j =
i; j < nfitpars; j++) {
6412 a(
i, j) - weightderiv(measno,
i) * weightderiv(measno, j)
6415 weightderiv(measno,
i) = 0;
6419 }
else if (trtrecal) {
6420 double *
errors = state->measurementErrors();
6421 double olderror =
errors[0];
6423 const auto *
const thisMeasurement{state->measurement()};
6429 if (oldrot->prepRawData() !=
nullptr) {
6430 double dcradius = oldrot->prepRawData()->localPosition()[
Trk::driftRadius];
6432 double trackradius = state->trackParameters()->parameters()[
Trk::driftRadius];
6434 std::unique_ptr<const Trk::RIO_OnTrack> newrot =
nullptr;
6435 double distance = std::abs(std::abs(trackradius) - dcradius);
6437 if (
distance < scalefactor * dcerror && (olderror > 1. || trackradius * oldradius < 0)) {
6438 newrot.reset(
m_ROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6439 }
else if (
distance > scalefactor * dcerror && olderror < 1.) {
6440 newrot.reset(
m_broadROTcreator->correct(*oldrot->prepRawData(), *state->trackParameters(), ctx));
6443 if (newrot !=
nullptr) {
6445 hitrecalibrated =
true;
6449 if ((measno < 0) or (measno >= (
int)
res.size())) {
6450 throw std::runtime_error(
6451 "'res' array index out of range in TrkGlobalChi2Fitter/src/GlobalChi2Fitter.cxx:" +
std::to_string(__LINE__)
6455 double oldres =
res[measno];
6456 double newres = newradius - state->trackParameters()->parameters()[
Trk::driftRadius];
6458 state->setMeasurement(std::move(newrot));
6460 trajectory.updateTRTHitCount(stateno, olderror);
6462 for (
int i = 0;
i < nfitpars;
i++) {
6463 if (weightderiv(measno,
i) == 0) {
6467 b[
i] -= weightderiv(measno,
i) * (oldres / olderror - (newres * olderror) / (newerror * newerror));
6469 for (
int j =
i; j < nfitpars; j++) {
6473 !cache.m_phiweight.empty() &&
6476 i < nperpars + 2 * nscats &&
6477 (
i - nperpars) % 2 == 0
6479 weight = cache.m_phiweight[(
i - nperpars) / 2];
6484 a(
i, j) + weightderiv(measno,
i) * weightderiv(measno, j) * ((olderror * olderror) / (newerror * newerror) - 1) *
weight
6487 weightderiv(measno,
i) *= olderror / newerror;
6490 res[measno] = newres;
6491 err[measno] = newerror;
6500 measno += state->numberOfMeasuredParameters();
6504 if (trajectory.nDOF() < 0) {
6509 if (outlierremoved || hitrecalibrated) {
6511 trajectory.setConverged(
false);
6513 cache.m_miniter =
it + 2;
◆ setMinIterations()
void Trk::GlobalChi2Fitter::setMinIterations |
( |
int |
| ) |
|
|
privatevirtual |
Definition at line 8443 of file GlobalChi2Fitter.cxx.
8445 (
"Configure the minimum number of Iterations via jobOptions");
◆ throwFailedToGetTrackingGeomtry()
void Trk::GlobalChi2Fitter::throwFailedToGetTrackingGeomtry |
( |
| ) |
const |
|
private |
◆ trackingGeometry()
Definition at line 1132 of file GlobalChi2Fitter.h.
1135 if (!cache.m_trackingGeometry)
1137 return cache.m_trackingGeometry;
◆ tryToConverge()
Definition at line 5408 of file GlobalChi2Fitter.cxx.
5415 const double oldChi2 = trajectory.prevchi2();
5416 const double newChi2 = trajectory.chi2();
5421 const double nDOF = trajectory.nDOF();
5422 const double oldRedChi2 = (nDOF > 0) ? oldChi2 / nDOF : 0;
5423 const double newRedChi2 = (nDOF > 0) ? newChi2 / nDOF : 0;
5426 trajectory.prefit() > 0 && (
5427 (newRedChi2 < 2 &&
it != 0) ||
5428 (newRedChi2 < oldRedChi2 + .1 && std::abs(newRedChi2 - oldRedChi2) < 1 &&
it != 1)
5431 trajectory.setConverged(
true);
5437 const int nsihits = trajectory.numberOfSiliconHits();
5438 const int ntrthits = trajectory.numberOfTRTHits();
5439 const int nhits = trajectory.numberOfHits();
5441 int miniter = (nsihits != 0 && nsihits + ntrthits == nhits) ? 1 : 2;
5442 miniter =
std::max(miniter, cache.m_miniter);
5444 if (
it >= miniter && std::abs(oldChi2 - newChi2) < 1) {
5445 trajectory.setConverged(
true);
◆ tryToWeightAfromMaterial()
Definition at line 5827 of file GlobalChi2Fitter.cxx.
5836 const int nPerPars = trajectory.numberOfPerigeeParameters();
5842 bool weightChanged =
false;
5849 double newPhiWeight = 1.1;
5850 double newThetaWeight = 1.001;
5851 if (trajectory.prefit() == 0) {
5857 newPhiWeight = 1.00000001;
5858 }
else if (
it == 1) {
5859 newPhiWeight = 1.0000001;
5860 }
else if (
it <= 3) {
5861 newPhiWeight = 1.0001;
5862 }
else if (
it <= 6) {
5863 newPhiWeight = 1.01;
5866 if (newRedChi2 > oldRedChi2 - 1 && newRedChi2 < oldRedChi2) {
5867 newPhiWeight = 1.0001;
5868 newThetaWeight = 1.0001;
5869 }
else if (newRedChi2 > oldRedChi2 - 25 && newRedChi2 < oldRedChi2) {
5870 newPhiWeight = 1.001;
5871 newThetaWeight = 1.0001;
5878 std::size_t scatno = 0;
5883 for (
const auto & state : trajectory.trackStates()) {
5884 const GXFMaterialEffects *meff = state->materialEffects();
5886 if (meff ==
nullptr) {
5890 const bool isValidPlaneSurface =
5892 static_cast<const PlaneSurface *
>(&state->associatedSurface()) !=
nullptr;
5897 if (meff->deltaE() == 0 || (trajectory.prefit() == 0 && isValidPlaneSurface)) {
5898 weightChanged =
true;
5900 const int scatNoIndex = 2 * scatno + nPerPars;
5902 if (trajectory.prefit() == 0 && meff->sigmaDeltaPhi() != 0) {
5903 if (scatno >= cache.m_phiweight.size()) {
5905 message <<
"scatno is out of range " << scatno <<
" !< " << cache.m_phiweight.size();
5906 throw std::range_error(
message.str());
5914 a(scatNoIndex, scatNoIndex) /= cache.m_phiweight[scatno];
5917 cache.m_phiweight[scatno] = newPhiWeight;
5918 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5919 }
else if (trajectory.prefit() >= 2) {
5920 a(scatNoIndex, scatNoIndex) *= newPhiWeight;
5921 a(scatNoIndex + 1, scatNoIndex + 1) *= newThetaWeight;
5933 meff->sigmaDeltaPhi() != 0 &&
5934 (trajectory.prefit() == 0 || meff->deltaE() == 0)
5948 trajectory.prefit() == 2 &&
5950 trajectory.numberOfBrems() > 0 &&
5951 (newRedChi2 < oldRedChi2 - 25 || newRedChi2 > oldRedChi2)
5956 return weightChanged;
◆ updateEnergyLoss()
Definition at line 7965 of file GlobalChi2Fitter.cxx.
7978 ATH_MSG_DEBUG(
"Angles out of range, phi: " << newphi <<
" theta: " << newtheta);
7982 double newqoverp = 0;
7984 if (meff.sigmaDeltaE() <= 0) {
7989 double newp2 = oldp * oldp -
sign * 2 * std::abs(meff.deltaE()) * std::sqrt(
mass *
mass + oldp * oldp) + meff.deltaE() * meff.deltaE();
7996 newqoverp = std::copysign(1 / std::sqrt(newp2),
old[
Trk::qOverP]);
8002 return surf.createUniqueTrackParameters(
8003 old[0],
old[1], newphi, newtheta, newqoverp, std::nullopt
◆ updateFitParameters()
Method to update peregee parameters, scattering angles, and brems.
Tries to solve the system [A] * deltaParameters = b and then update in the trajectory all parameters used for the fit. Returns also a status.
Definition at line 6132 of file GlobalChi2Fitter.cxx.
6145 Eigen::LLT<Eigen::MatrixXd> llt(lu_m);
6147 if (llt.info() != Eigen::Success) {
6157 const int nscat = trajectory.numberOfScatterers();
6158 const int nbrem = trajectory.numberOfBrems();
6159 const int nperparams = trajectory.numberOfPerigeeParameters();
6171 double d0 = refpar->parameters()[
Trk::d0];
6172 double z0 = refpar->parameters()[
Trk::z0];
6175 double qoverp = refpar->parameters()[
Trk::qOverP];
6177 if (nperparams > 0) {
6178 d0 += deltaParameters[0];
6179 z0 += deltaParameters[1];
6180 phi += deltaParameters[2];
6181 theta += deltaParameters[3];
6182 qoverp = (trajectory.m_straightline) ? 0 : .001 * deltaParameters[4] + qoverp;
6194 std::vector < std::pair < double, double >>&scatangles = trajectory.scatteringAngles();
6195 for (
int i = 0;
i < nscat;
i++) {
6196 scatangles[
i].first += deltaParameters[2 *
i + nperparams];
6197 scatangles[
i].second += deltaParameters[2 *
i + nperparams + 1];
6203 std::vector < double >&delta_ps = trajectory.brems();
6204 for (
int i = 0;
i < nbrem;
i++) {
6205 delta_ps[
i] += deltaParameters[nperparams + 2 * nscat +
i];
6211 std::unique_ptr<const TrackParameters> newper(
6212 trajectory.referenceParameters()->associatedSurface().createUniqueTrackParameters(
6220 trajectory.setReferenceParameters(std::move(newper));
6221 trajectory.setScatteringAngles(scatangles);
6222 trajectory.setBrems(delta_ps);
◆ updatePixelROTs()
Update the Pixel ROT using the current trajectory/local track parameters.
Definition at line 6227 of file GlobalChi2Fitter.cxx.
6233 if ( trajectory.numberOfSiliconHits() == 0) {
6242 if (!splitProbContainer.isValid()) {
6246 std::vector<std::unique_ptr<GXFTrackState>> &
states = trajectory.trackStates();
6249 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
6250 int nfitpars = trajectory.numberOfFitParameters();
6253 for (
size_t stateno = 0; stateno <
states.size(); stateno++) {
6258 measno +=
states[stateno-1]->numberOfMeasuredParameters();
6261 std::unique_ptr<GXFTrackState> & state =
states[stateno];
6271 const PrepRawData *prd{};
6273 const auto *
const rot =
static_cast<const RIO_OnTrack *
>(pMeas);
6274 prd = rot->prepRawData();
6284 const auto &splitProb = splitProbContainer->splitProbability(
pixelCluster);
6285 if (!splitProb.isSplit()) {
6286 ATH_MSG_DEBUG(
"Pixel cluster is not split so no need to update" );
6290 std::unique_ptr < const RIO_OnTrack > newrot;
6291 double *olderror = state->measurementErrors();
6294 double newerror[5] = {-1,-1,-1,-1,-1};
6295 double newres[2] = {-1,-1};
6297 newrot.reset(
m_ROTcreator->correct(*prd, *trackpars, evtctx));
6304 newerror[0] = std::sqrt(
covmat(0, 0));
6305 newres[0] = newrot->localParameters()[
Trk::locX] - trackpars->parameters()[
Trk::locX];
6306 newerror[1] = std::sqrt(
covmat(1, 1));
6307 newres[1] = newrot->localParameters()[
Trk::locY] - trackpars->parameters()[
Trk::locY];
6309 if (
a.cols() != nfitpars) {
6314 for(
int k =0;
k<2;
k++ ){
6315 double oldres =
res[measno+
k];
6316 res[measno+
k] = newres[
k];
6317 err[measno+
k] = newerror[
k];
6319 for (
int i = 0;
i < nfitpars;
i++) {
6320 if (weightderiv(measno+
k,
i) == 0) {
6324 b[
i] -= weightderiv(measno+
k,
i) * (oldres / olderror[
k] - (newres[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]));
6326 for (
int j =
i; j < nfitpars; j++) {
6330 weightderiv(measno+
k,
i) *
6331 weightderiv(measno+
k, j) *
6332 ((olderror[
k] * olderror[
k]) / (newerror[
k] * newerror[
k]) - 1)
6336 weightderiv(measno+
k,
i) *= olderror[
k] / newerror[
k];
6340 state->setMeasurement(std::move(newrot));
6341 state->setMeasurementErrors(newerror);
◆ updateSystemWithMaxBremPull()
Definition at line 5449 of file GlobalChi2Fitter.cxx.
5457 if (state_maxbrempull ==
nullptr) {
5461 state_maxbrempull->materialEffects()->setSigmaDeltaE(
5462 10 * state_maxbrempull->materialEffects()->sigmaDeltaEPos()
5465 state_maxbrempull->materialEffects()->setKink(
true);
5467 const int nbrem = trajectory.numberOfBrems();
5469 const int nmeas = (
int)
res.size();
5472 const double oldError =
error[nmeas - nbrem + bremno_maxbrempull];
5473 const double newError = .001 * state_maxbrempull->materialEffects()->sigmaDeltaE();
5474 error[nmeas - nbrem + bremno_maxbrempull] = newError;
5476 const int nFitPars = trajectory.numberOfFitParameters();
5477 if (
a.cols() != nFitPars) {
5481 const double errorRatio = oldError / newError;
5482 const double errorReductionRatio = 1 -
std::pow(errorRatio, 2);
5484 Amg::MatrixX & weightderiv = trajectory.weightedResidualDerivatives();
5485 for (
int i = 0;
i < nFitPars;
i++) {
5486 if (weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) == 0) {
5490 for (
int j =
i; j < nFitPars; j++) {
5491 const double newaij =
a(
i, j) - errorReductionRatio *
5492 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *
5493 weightderiv(nmeas - nbrem + bremno_maxbrempull, j);
5495 a.fillSymmetric(
i, j, newaij);
5497 weightderiv(nmeas - nbrem + bremno_maxbrempull,
i) *= errorRatio;
◆ ATLAS_THREAD_SAFE
std::array<std::atomic<unsigned int>, S_MAX_VALUE> m_fit_status Trk::GlobalChi2Fitter::ATLAS_THREAD_SAFE = {} |
|
mutableprivate |
◆ m_acceleration
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_acceleration {this, "Acceleration", false} |
|
private |
◆ m_asymeloss
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_asymeloss {this, "AsymmetricEnergyLoss", true} |
|
private |
◆ m_boundaryCheckTool
ToolHandle<IBoundaryCheckTool> Trk::GlobalChi2Fitter::m_boundaryCheckTool {this, "BoundaryCheckTool", "", "Boundary checking tool for detector sensitivities" } |
|
private |
◆ m_broadROTcreator
ToolHandle<IRIO_OnTrackCreator> Trk::GlobalChi2Fitter::m_broadROTcreator {this, "BroadRotCreatorTool", "", ""} |
|
private |
◆ m_calomat
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_calomat {this, "MuidMat", false} |
|
private |
◆ m_caloMaterialProvider
◆ m_calotool
◆ m_calotoolparam
◆ m_chi2cut
Gaudi::Property<double> Trk::GlobalChi2Fitter::m_chi2cut {this, "TrackChi2PerNDFCut", 1.e15} |
|
private |
◆ m_clusterSplitProbContainer
◆ m_createSummary
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_createSummary {this, "CreateTrackSummary", true} |
|
private |
◆ m_decomposesegments
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_decomposesegments {this, "DecomposeSegments", true} |
|
private |
◆ m_DetID
◆ m_domeastrackpar
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_domeastrackpar {this, "MeasuredTrackParameters", true} |
|
private |
◆ m_elosstool
◆ m_extensioncuts
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_extensioncuts {this, "TRTExtensionCuts", true} |
|
private |
◆ m_extmat
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_extmat {this, "ExtrapolatorMaterial", true} |
|
private |
◆ m_extrapolator
◆ m_field_cache_key
Initial value:{
this,
"AtlasFieldCacheCondObj",
"fieldCondObj",
"Trk::GlobalChi2Fitter field conditions object key"
}
Definition at line 1157 of file GlobalChi2Fitter.h.
◆ m_fillderivmatrix
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_fillderivmatrix {this, "FillDerivativeMatrix", false} |
|
private |
◆ m_fiteloss
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_fiteloss {this, "FitEnergyLoss", false} |
|
private |
◆ m_fixbrem
Gaudi::Property<int> Trk::GlobalChi2Fitter::m_fixbrem {this, "FixBrem", -1} |
|
private |
◆ m_getmaterialfromtrack
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_getmaterialfromtrack {this, "GetMaterialFromTrack", true} |
|
private |
◆ m_holeSearch
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_holeSearch {this, "DoHoleSearch", false} |
|
private |
◆ m_idVolume
◆ m_kinkfinding
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_kinkfinding {this, "KinkFinding", false} |
|
private |
◆ m_matupdator
◆ m_maxit
Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxit {this, "MaxIterations", 30} |
|
private |
◆ m_maxitPixelROT
Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxitPixelROT {this, "IterationsToRebuildPixelRots", 0} |
|
private |
◆ m_maxoutliers
Gaudi::Property<int> Trk::GlobalChi2Fitter::m_maxoutliers {this, "MaxOutliers", 10} |
|
private |
◆ m_miniter
Gaudi::Property<int> Trk::GlobalChi2Fitter::m_miniter {this, "MinimumIterations", 1} |
|
private |
◆ m_minphfcut
Gaudi::Property<double> Trk::GlobalChi2Fitter::m_minphfcut {this, "MinPHFCut", 0.} |
|
private |
◆ m_navigator
ToolHandle<INavigator> Trk::GlobalChi2Fitter::m_navigator {this, "NavigatorTool", "Trk::Navigator/CosmicsNavigator", ""} |
|
private |
◆ m_numderiv
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_numderiv {this, "NumericalDerivs", false} |
|
private |
◆ m_outlcut
Gaudi::Property<double> Trk::GlobalChi2Fitter::m_outlcut {this, "OutlierCut", 5.0} |
|
private |
◆ m_p
Gaudi::Property<double> Trk::GlobalChi2Fitter::m_p {this, "Momentum", 0.0} |
|
private |
◆ m_propagator
ToolHandle<IPropagator> Trk::GlobalChi2Fitter::m_propagator {this, "PropagatorTool", "", ""} |
|
private |
◆ m_redoderivs
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_redoderivs {this, "RecalculateDerivatives", false} |
|
private |
◆ m_reintoutl
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_reintoutl {this, "ReintegrateOutliers", false} |
|
private |
◆ m_rejectLargeNScat
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_rejectLargeNScat {this, "RejectLargeNScat", false} |
|
private |
◆ m_residualPullCalculator
◆ m_ROTcreator
ToolHandle<IRIO_OnTrackCreator> Trk::GlobalChi2Fitter::m_ROTcreator {this, "RotCreatorTool", "", ""} |
|
private |
◆ m_scalefactor
Gaudi::Property<double> Trk::GlobalChi2Fitter::m_scalefactor {this, "TRTTubeHitCut", 2.5} |
|
private |
◆ m_scattool
◆ m_signedradius
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_signedradius {this, "SignedDriftRadius", true} |
|
private |
◆ m_sirecal
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_sirecal {this, "RecalibrateSilicon", false} |
|
private |
◆ m_storemat
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_storemat {this, "StoreMaterialOnTrack", true} |
|
private |
◆ m_straightlineprop
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_straightlineprop {this, "StraightLine", true} |
|
private |
◆ m_trackingGeometryReadKey
Initial value:{
this,
"TrackingGeometryReadKey",
"AtlasTrackingGeometry",
"Key of the TrackingGeometry conditions data."
}
Definition at line 1150 of file GlobalChi2Fitter.h.
◆ m_trtrecal
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_trtrecal {this, "RecalibrateTRT", false} |
|
private |
◆ m_updator
ToolHandle<IUpdator> Trk::GlobalChi2Fitter::m_updator {this, "MeasurementUpdateTool", "", ""} |
|
private |
◆ m_useCaloTG
Gaudi::Property<bool> Trk::GlobalChi2Fitter::m_useCaloTG {this, "UseCaloTG", false} |
|
private |
The documentation for this class was generated from the following files:
def retrieve(aClass, aKey=None)
bool is_pixel(Identifier id) const
std::unique_ptr< const TrackParameters > makePerigee(Cache &, const TrackParameters &, const ParticleHypothesis) const
void makeProtoState(Cache &, GXFTrajectory &, const TrackStateOnSurface *, int index=-1) const
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const TrackingGeometry * retrieveTrackingGeometry(const EventContext &ctx) const
bool is_rpc(Identifier id) const
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Gaudi::Property< bool > m_createSummary
Const iterator class for DataVector/DataList.
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
ToolHandle< INavigator > m_navigator
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
void fillResidualsAndErrors(const EventContext &ctx, const Cache &cache, GXFTrajectory &trajectory, const int it, Amg::VectorX &b, int &bremno_maxbrempull, GXFTrackState *&state_maxbrempull) const
bool tryToWeightAfromMaterial(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a, const bool doDeriv, const int it, const double oldRedChi2, const double newRedChi2) const
@ DeadElement
outside the element
Gaudi::Property< double > m_outlcut
@ z
global position (cartesian)
bool ensureValidEntranceMuonSpectrometer(const EventContext &ctx, Cache &cache) const
Gaudi::Property< bool > m_numderiv
ToolHandle< IMaterialEffectsOnTrackProvider > m_calotool
std::string find(const std::string &s)
return a remapped string
ToolHandle< IMaterialEffectsUpdator > m_matupdator
bool is_csc(Identifier id) const
Scalar perp() const
perp method - perpenticular length
bool inside(const Amg::Vector3D &gp, double tol=0.) const
Inside() method for checks.
Gaudi::Property< bool > m_holeSearch
bool is_sct(Identifier id) const
const Amg::Vector3D & position() const
Access method for the position.
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
FitterStatusCode runIteration(const EventContext &ctx, Cache &cache, GXFTrajectory &trajectory, const int it, Amg::SymMatrixX &a, Amg::VectorX &b, Amg::SymMatrixX &lu, bool &doDeriv) const
@ numberOfSCTDeadSensors
number of TRT hits
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
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
void fillFirstLastMeasurement(Cache &cache, GXFTrajectory &trajectory) const
represents a deflection of the track caused through multiple scattering in material.
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
ToolHandle< IPropagator > m_propagator
@ loc2
generic first and second local coordinate
void runTrackCleanerTRT(Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, Amg::SymMatrixX &, bool, bool, int, const EventContext &ctx) const
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
void addMaterialUpdateTrajectory(Cache &cache, GXFTrajectory &track, int offset, std::vector< std::pair< const Layer *, const Layer * >> &layers, const TrackParameters *ref1, const TrackParameters *ref2, ParticleHypothesis mat) const
Given layer information, probe those layers for scatterers and add them to a track.
std::vector< size_t > vec
bool is_trt(Identifier id) const
std::vector< std::unique_ptr< TrackParameters > > holesearchExtrapolation(const EventContext &ctx, const TrackParameters &src, const GXFTrackState &dst, PropDirection propdir) const
Helper method which performs an extrapolation with additional logic for hole search.
#define ATH_MSG_VERBOSE(x)
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 ...
void updateSystemWithMaxBremPull(GXFTrajectory &trajectory, const int bremno_maxbrempull, GXFTrackState *state_maxbrempull, Amg::SymMatrixX &a) const
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
const LayerArray * confinedLayers() const
Return the subLayer array.
Gaudi::Property< int > m_maxit
@ StraightTrack
A straight track.
AmgSymMatrix(5) &GXFTrackState
@ SlimmedTrack
A slimmed track.
represents the full description of deflection and e-loss of a track in material.
ToolHandle< IResidualPullCalculator > m_residualPullCalculator
std::vector< std::reference_wrapper< GXFTrackState > > holeSearchStates(GXFTrajectory &trajectory) const
Extracts a collection of track states which are important for hole search.
static void addMaterialGetLayers(Cache &cache, std::vector< std::pair< const Layer *, const Layer * >> &layers, std::vector< std::pair< const Layer *, const Layer * >> &uplayers, const std::vector< std::unique_ptr< GXFTrackState >> &states, GXFTrackState &first, GXFTrackState &last, const TrackParameters *refpar, bool hasmat)
Collect all possible layers that a given track could have passed through.
static void makeTrackFillDerivativeMatrix(Cache &, GXFTrajectory &)
ToolHandle< Trk::ITrkMaterialProviderTool > m_caloMaterialProvider
double getDistance(const xAOD::Vertex *vtx1, const xAOD::Vertex *vtx2)
void compensatePhiWeights(Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) 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
void tryToConverge(const Cache &cache, GXFTrajectory &trajectory, const int it) const
ToolHandle< IMultipleScatteringUpdator > m_scattool
std::unique_ptr< Track > makeTrack(const EventContext &ctx, Cache &, GXFTrajectory &, const ParticleHypothesis) const
@ FullField
Field is set to be realistic, but within a given Volume.
bool processTrkVolume(Cache &, const Trk::TrackingVolume *tvol) const
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
@ BremPoint
This represents a brem point on the track, and so will contain TrackParameters and MaterialEffectsBas...
std::optional< TransportJacobian > numericalDerivatives(const EventContext &ctx, const TrackParameters *, const Surface &, PropDirection, const MagneticFieldProperties &) const
@ numberOfPixelHoles
number of pixels which have a ganged ambiguity.
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
Gaudi::Property< bool > m_acceleration
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
@ NoField
Field is set to 0., 0., 0.,.
Ensure that the ATLAS eigen extensions are properly loaded.
void addMaterial(const EventContext &ctx, Cache &, GXFTrajectory &, const TrackParameters *, ParticleHypothesis) const
void updatePixelROTs(GXFTrajectory &, Amg::SymMatrixX &, Amg::VectorX &, const EventContext &evtctx) const
Update the Pixel ROT using the current trajectory/local track parameters.
std::vector< std::string > remainder(const std::vector< std::string > &v1, const std::vector< std::string > &v2)
std::string to_string(const DetectorType &type)
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
GXFTrajectory * runTrackCleanerSilicon(const EventContext &ctx, Cache &, GXFTrajectory &, Amg::SymMatrixX &, Amg::SymMatrixX &, Amg::VectorX &, bool) const
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)
void fillBfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::VectorX &b) const
Eigen::Matrix< double, 3, 1 > Vector3D
FitterStatusCode updateFitParameters(GXFTrajectory &, const Amg::VectorX &, const Amg::SymMatrixX &) const
Method to update peregee parameters, scattering angles, and brems.
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
void fillAfromMeasurements(const Cache &cache, GXFTrajectory &trajectory, Amg::SymMatrixX &a) const
const Amg::Vector3D & momentum() const
Access method for the momentum.
Track * backupCombinationStrategy(const EventContext &ctx, Cache &, const Track &, const Track &, GXFTrajectory &, std::vector< MaterialEffectsOnTrack > &) const
@ Biased
RP with track state including the hit.
bool is_muon(Identifier id) const
@ GlobalChi2Fitter
Track's from Thijs' global chi^2 fitter.
void makeProtoStateFromMeasurement(Cache &, GXFTrajectory &, const MeasurementBase *, const TrackParameters *trackpar=nullptr, bool isoutlier=false, int index=-1) const
Gaudi::Property< double > m_p
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
#define ATH_MSG_WARNING(x)
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
def delta_phi(phi1, phi2)
Gaudi::Property< int > m_maxitPixelROT
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
Track * fitIm(const EventContext &ctx, Cache &cache, const Track &inputTrack, const RunOutlierRemoval runOutlier, const ParticleHypothesis matEffects) const
std::unique_ptr< T > unique_clone(const T *v)
Gaudi::Property< double > m_scalefactor
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
ToolHandle< IUpdator > m_updator
void addIDMaterialFast(const EventContext &ctx, Cache &cache, GXFTrajectory &track, const TrackParameters *parameters, ParticleHypothesis part) const
A faster strategy for adding scatter material to tracks, works only for inner detector tracks.
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const PrepRawDataSet &, const TrackParameters &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=nonInteracting) const override final
@ PseudoMeasurementOnTrack
bool is_stgc(Identifier id) const
Gaudi::Property< bool > m_calomat
Gaudi::Property< bool > m_rejectLargeNScat
Gaudi::Property< bool > m_trtrecal
bool consistentSurfaces(U)
virtual MaterialEffectsDerivedType derivedType() const =0
Returns the concrete derived type.
constexpr int pow(int base, int exp) noexcept
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
bool ensureValidEntranceCalo(const EventContext &ctx, Cache &cache) const
Scalar mag() const
mag method
std::variant< std::unique_ptr< const TrackParameters >, FitterStatusCode > updateEnergyLoss(const Surface &, const GXFMaterialEffects &, const TrackParameters &, double, int) const
void fillAfromScatterers(GXFTrajectory &trajectory, Amg::SymMatrixX &a) const
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
static void calculateDerivatives(GXFTrajectory &)
float nSiHits(const U &p)
bool is_mdt(Identifier id) const
Gaudi::Property< bool > m_useCaloTG
void fillDerivatives(GXFTrajectory &traj) const
virtual double r() const override final
This method returns the radius.
@ OutlierLogicFailure
outlier logic failed
SG::ReadHandleKey< Trk::ClusterSplitProbabilityContainer > m_clusterSplitProbContainer
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > SymMatrixX
@ Unknown
Track fitter not defined.