From 9755ce18a00ea3f562c749023328b6e1a56cf9e1 Mon Sep 17 00:00:00 2001 From: shahoian Date: Sun, 14 Jun 2026 14:36:59 +0200 Subject: [PATCH 1/2] Changes in TPC residuals validation The meaning of the scdcalib.writeUnfiltered (default: false) ConfigurableParam has changed. If set to true, tracks that do not pass the validateTrack check will no longer be discarded but will instead be stored in the standard unbinned residuals data. The int8_t TrackData::filterFlag field provides the status of the validation check: -1: no validation check was requested; 0: validation passed; >0: reason for failing validation (see the validateTrack method). A new ConfigurableParam, scdcalib.writeValidationData (default: false), has been added. If set to true, the TrackValidationData object, together with the corresponding TrackData, will be written to a dedicated debug tree, "valdata", in "track_interpolation_dbg.root" (or "track_interpolation_dbg_.root" when multiple lanes are used). If a track passes validation but some residuals are flagged as bad, the new ConfigurableParam scdcalib.keepRejectedResiduals (default: false) allows these residuals to be written to the unbinned residuals tree with the bool UnbinnedResid::rejected flag set (previously, such residuals were discarded). --- .../src/TPCInterpolationSpec.cxx | 18 +- .../src/TPCResidualWriterSpec.cxx | 2 - .../SpacePoints/SpacePointsCalibConfParam.h | 4 +- .../include/SpacePoints/TrackInterpolation.h | 111 +++-- .../include/SpacePoints/TrackResiduals.h | 3 +- .../SpacePoints/src/SpacePointCalibLinkDef.h | 3 + .../SpacePoints/src/TrackInterpolation.cxx | 422 ++++++++++-------- .../SpacePoints/src/TrackResiduals.cxx | 101 ++++- 8 files changed, 403 insertions(+), 261 deletions(-) diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx index 4f17aec3d49d5..fe2677416edc2 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx @@ -32,6 +32,7 @@ #include "SpacePoints/SpacePointsCalibConfParam.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/ControlService.h" +#include "Framework/DeviceSpec.h" using namespace o2::framework; using namespace o2::globaltracking; @@ -55,6 +56,9 @@ void TPCInterpolationDPL::init(InitContext& ic) if (mProcessSeeds && mSources != mSourcesMap) { LOG(fatal) << "process-seeds option is not compatible with using different track sources for vDrift and map extraction"; } + int lane = ic.services().get().inputTimesliceId; + int maxLanes = ic.services().get().maxInputTimeslices; + mInterpolation.setLane(lane, maxLanes); } void TPCInterpolationDPL::updateTimeDependentParams(ProcessingContext& pc) @@ -136,13 +140,6 @@ void TPCInterpolationDPL::run(ProcessingContext& pc) mInterpolation.process(); mTimer.Stop(); LOGF(info, "TPC interpolation timing: Cpu: %.3e Real: %.3e s", mTimer.CpuTime(), mTimer.RealTime()); - if (SpacePointsCalibConfParam::Instance().writeUnfiltered) { - // these are the residuals and tracks before outlier rejection; they are not used in production - pc.outputs().snapshot(Output{"GLO", "TPCINT_RES", 0}, mInterpolation.getClusterResidualsUnfiltered()); - if (mSendTrackData) { - pc.outputs().snapshot(Output{"GLO", "TPCINT_TRK", 0}, mInterpolation.getReferenceTracksUnfiltered()); - } - } pc.outputs().snapshot(Output{"GLO", "UNBINNEDRES", 0}, mInterpolation.getClusterResiduals()); pc.outputs().snapshot(Output{"GLO", "DETINFORES", 0}, mInterpolation.getClusterResidualsDetInfo()); pc.outputs().snapshot(Output{"GLO", "TRKREFS", 0}, mInterpolation.getTrackDataCompact()); @@ -157,6 +154,7 @@ void TPCInterpolationDPL::run(ProcessingContext& pc) void TPCInterpolationDPL::endOfStream(EndOfStreamContext& ec) { + mInterpolation.finalize(); LOGF(info, "TPC residuals extraction total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } @@ -183,12 +181,6 @@ DataProcessorSpec getTPCInterpolationSpec(GTrackID::mask_t srcCls, GTrackID::mas dataRequest->inputs, true); o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs); - if (SpacePointsCalibConfParam::Instance().writeUnfiltered) { - outputs.emplace_back("GLO", "TPCINT_TRK", 0, Lifetime::Timeframe); - if (sendTrackData) { - outputs.emplace_back("GLO", "TPCINT_RES", 0, Lifetime::Timeframe); - } - } outputs.emplace_back("GLO", "UNBINNEDRES", 0, Lifetime::Timeframe); outputs.emplace_back("GLO", "DETINFORES", 0, Lifetime::Timeframe); outputs.emplace_back("GLO", "TRKREFS", 0, Lifetime::Timeframe); diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx index 8b06444bdb9b3..68c9b3abfe4ef 100644 --- a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCResidualWriterSpec.cxx @@ -35,8 +35,6 @@ DataProcessorSpec getTPCResidualWriterSpec(bool writeTrackData, bool debugOutput return MakeRootTreeWriterSpec("tpc-residuals-writer", "o2residuals_tpc.root", "residualsTPC", - BranchDefinition>{InputSpec{"tracksUnfiltered", "GLO", "TPCINT_TRK", 0}, "tracksUnfiltered", ((writeUnfiltered && writeTrackData) ? 1 : 0)}, - BranchDefinition>{InputSpec{"residualsUnfiltered", "GLO", "TPCINT_RES", 0}, "residualsUnfiltered", (writeUnfiltered ? 1 : 0)}, BranchDefinition>{InputSpec{"residuals", "GLO", "UNBINNEDRES"}, "residuals"}, BranchDefinition>{InputSpec{"detInfo", "GLO", "DETINFORES"}, "detInfo"}, BranchDefinition>{InputSpec{"trackRefs", "GLO", "TRKREFS"}, "trackRefs"}, diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h index 8b884209dd697..5e25b4717118c 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h @@ -63,7 +63,8 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper0 giving the reason of rejection + bool keepRejectedResiduals{false}; ///< if set, keep rejected residuals setting rejected flag int nMALong{15}; ///< number of points to be used for moving average (long range) int nMAShort{3}; ///< number of points to be used for estimation of distance from local line (short range) float maxRejFrac{.15f}; ///< if the fraction of rejected clusters of a track is higher, the full track is invalidated @@ -74,6 +75,7 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper(dyIn * 0x7fff / param::MaxResid)), - dz(static_cast(dzIn * 0x7fff / param::MaxResid)), - tgSlp(static_cast(tgSlpIn * 0x7fff / param::MaxTgSlp)), - y(static_cast(yIn * 0x7fff / param::MaxY)), - z(static_cast(zIn * 0x7fff / param::MaxZ)), - row(rowIn), - sec(secIn), - channel(chanIn) {} + UnbinnedResid(float dyIn, float dzIn, float tgSlpIn, float yIn, float zIn, + unsigned char rowIn, unsigned char secIn, short chanIn = -1, bool rejFlag = false) : dy(static_cast(dyIn * 0x7fff / param::MaxResid)), + dz(static_cast(dzIn * 0x7fff / param::MaxResid)), + tgSlp(static_cast(tgSlpIn * 0x7fff / param::MaxTgSlp)), + y(static_cast(yIn * 0x7fff / param::MaxY)), + z(static_cast(zIn * 0x7fff / param::MaxZ)), + row(rowIn), + sec(secIn), + channel(chanIn), + rejected(rejFlag) {} short dy{0}; ///< residual in y short dz{0}; ///< residual in z short tgSlp{0}; ///< tan of the phi angle between padrow and track @@ -88,6 +91,7 @@ struct UnbinnedResid { unsigned char row{0}; ///< TPC pad row unsigned char sec{0}; ///< TPC sector (0..35) short channel{-1}; ///< extra channel info (ITS chip ID, TRD chamber, TOF main pad within the sector) + bool rejected{false}; ///< residual is flagged as rejected in the validateTrack bool isTPC() const { return row < constants::MAXGLOBALPADROW; } bool isTRD() const { return row >= 160 && row < 166; } @@ -103,7 +107,7 @@ struct UnbinnedResid { static void checkInitDone(); static bool gInitDone; - ClassDefNV(UnbinnedResid, 2); + ClassDefNV(UnbinnedResid, 3); }; struct DetInfoResid { // detector info associated with residual @@ -153,12 +157,13 @@ struct DetInfoResid { // detector info associated with residual /// Structure for the information required to associate each residual with a given track type (ITS-TPC-TRD-TOF, etc) struct TrackDataCompact { TrackDataCompact() = default; - TrackDataCompact(uint32_t idx, std::array mlt, uint8_t nRes, uint8_t source, uint8_t nextraRes = 0) : idxFirstResidual(idx), multStack{mlt}, nResiduals(nRes), sourceId(source), nExtDetResid(nextraRes) {} + TrackDataCompact(uint32_t idx, std::array mlt, uint8_t nRes, uint8_t source, uint8_t nextraRes = 0, int8_t filt = -1) : idxFirstResidual(idx), multStack{mlt}, nResiduals(nRes), sourceId(source), filterFlag(filt), nExtDetResid(nextraRes) {} uint32_t idxFirstResidual; ///< the index of the first residual from this track std::array multStack{}; // multiplicity in the stack packed as asinh(x*0.05)/0.05 uint8_t nResiduals; ///< total number of TPC residuals associated to this track uint8_t nExtDetResid = 0; ///< number of external detectors (wrt TPC) residuals stored, on top of clIdx.getEntries uint8_t sourceId; ///< source ID obtained from the global track ID + int8_t filterFlag = -1; ///< -1: validation was not done, 0: validated, >0 : reason not passing validation, see validateTrack method void setMultStack(float v, int stack) { @@ -171,7 +176,7 @@ struct TrackDataCompact { } float getMultStackPacked(int stack) const { return multStack[stack]; } - ClassDefNV(TrackDataCompact, 3); + ClassDefNV(TrackDataCompact, 4); }; // TODO add to UnbinnedResid::sec flag if cluster was used or not @@ -191,7 +196,8 @@ struct TrackDataExtended { o2::tof::Cluster clsTOF{}; ///< the TOF cluster (if available) o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track uint8_t nExtDetResid = 0; ///< number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries - ClassDefNV(TrackDataExtended, 3); + int8_t filterFlag = -1; ///< -1: validation was not done, 0: validated, >0 : reason not passing validation, see validateTrack method + ClassDefNV(TrackDataExtended, 4); }; /// Structure filled for each track with track quality information and a vector with TPCClusterResiduals @@ -210,6 +216,7 @@ struct TrackData { unsigned short clAvailTOF{}; ///< whether or not track seed has a matched TOF cluster, if so, gives the resolution of the T0 in ps short TRDTrkltSlope[6] = {}; ///< TRD tracklet slope 0x7fff / param::MaxTRDSlope uint8_t nExtDetResid = 0; ///< number of external detectors (to TPC) residuals stored, on top of clIdx.getEntries + int8_t filterFlag = -1; ///< -1: validation was not done, 0: validated, >0 : reason not passing validation, see validateTrack method o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of TPC cluster residuals of this track std::array multStack{}; // multiplicity in the stack packed as asinh(x*0.05)/0.05 float getT0Error() const { return float(clAvailTOF); } @@ -226,7 +233,7 @@ struct TrackData { } float getMultStackPacked(int stack) const { return multStack[stack]; } - ClassDefNV(TrackData, 11); + ClassDefNV(TrackData, 12); }; /// \class TrackInterpolation @@ -245,6 +252,8 @@ class TrackInterpolation TrackInterpolation(const TrackInterpolation&) = delete; TrackInterpolation& operator=(const TrackInterpolation&) = delete; + ~TrackInterpolation(); + /// Enumeration for indexing the arrays of the CacheStruct enum { ExtOut = 0, ///< extrapolation outwards of ITS track @@ -269,16 +278,40 @@ class TrackInterpolation }; /// Structure for on-the-fly re-calculated track parameters at the validation stage - struct TrackParams { - TrackParams() = default; + struct ValidationPoint { + float xTrk{0.f}; + float yTrk{0.f}; + float zTrk{0.f}; + float xLab{0.f}; + float yLab{0.f}; + float sPath{0.f}; + float dy{0.f}; + float dz{0.f}; + float tglArr{0.f}; + float residHelixY{0.f}; + float residHelixZ{0.f}; + float diffYSmooth{0.f}; + float diffZSmooth{0.f}; + int8_t sec{0}; + bool flagRej{false}; + ClassDefNV(ValidationPoint, 1); + }; + + struct TrackValidationData { float qpt{0.f}; float tgl{0.f}; - std::array zTrk{}; - std::array xTrk{}; - std::array dy{}; - std::array dz{}; - std::array tglArr{}; - std::bitset flagRej{}; + float xcLab{0.f}; + float ycLab{0.f}; + float r{0.f}; + float zOffs{0.f}; + uint8_t nRej = 0; + std::vector points; + void clear() + { + points.clear(); + nRej = 0; + } + ClassDefNV(TrackValidationData, 1); }; // -------------------------------------- processing functions -------------------------------------------------- @@ -319,26 +352,26 @@ class TrackInterpolation /// Validates the given input track and its residuals /// \param trk The track parameters, e.g. q/pT, eta, ... /// \param params Structure with per pad information recalculated on the fly - /// \return true if the track could be validated, false otherwise - bool validateTrack(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const; + /// \return 0 if the track could be validated, otherwise returns rejection code + int8_t validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes); /// Filter out individual outliers from all cluster residuals of given track /// \return true for tracks which pass the cuts on e.g. max. masked clusters and false for rejected tracks - bool outlierFiltering(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const; + bool outlierFiltering(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes); /// Is called from outlierFiltering() and does the actual calculations (moving average filter etc.) /// \return The RMS of the long range moving average - float checkResiduals(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const; + float checkResiduals(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes); /// Calculates the differences in Y and Z for a given set of clusters to a fitted helix. /// First a circular fit in the azimuthal plane is performed and subsequently a linear fit in the transversal plane - bool compareToHelix(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const; + bool compareToHelix(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes); /// For a given set of points, calculate the differences from each point to the fitted lines from all other points in their neighbourhoods (+- nMAShort points) - void diffToLocLine(const int np, int idxOffset, const std::array& x, const std::array& y, std::array& diffY) const; + void diffToLocLine(TrackValidationData& params, int start, int np); /// For a given set of points, calculate their deviation from the moving average (build from the neighbourhood +- nMALong points) - void diffToMA(const int np, const std::array& y, std::array& diffMA) const; + void diffToMA(const int np, const std::array& y, std::array& diffMA); // -------------------------------------- settings -------------------------------------------------- void setNHBPerTF(int n) { mNHBPerTF = n; } @@ -381,8 +414,14 @@ class TrackInterpolation std::vector& getTrackDataCompact() { return mTrackDataCompact; } std::vector& getTrackDataExtended() { return mTrackDataExtended; } std::vector& getReferenceTracks() { return mTrackData; } - std::vector& getClusterResidualsUnfiltered() { return mClResUnfiltered; } - std::vector& getReferenceTracksUnfiltered() { return mTrackDataUnfiltered; } + + void setLane(int lID, int nL) + { + mLaneID = lID; + mNLanes = nL; + } + + void finalize(); private: static constexpr float sFloatEps{1.e-7f}; ///< float epsilon for robust linear fitting @@ -391,6 +430,8 @@ class TrackInterpolation // parameters + settings const SpacePointsCalibConfParam* mParams = nullptr; std::shared_ptr mTPCParam = nullptr; + int mLaneID = 0; + int mNLanes = 1; int mNHBPerTF = 32; int mNTPCOccBinLength = 16; ///< TPC occupancy bin length in TB float mNTPCOccBinLengthInv = 1.f / 16; ///< its inverse @@ -433,13 +474,14 @@ class TrackInterpolation std::vector mTrackDataExtended{}; ///< full tracking information for debugging std::vector mClRes{}; ///< residuals for each available TPC cluster of all tracks std::vector mDetInfoRes{}; ///< packed detector info associated with each residual - std::vector mTrackDataUnfiltered{}; ///< same as mTrackData, but for all tracks before outlier filtering - std::vector mClResUnfiltered{}; ///< same as mClRes, but for all residuals before outlier filtering + std::unique_ptr mDBGOut; // cache std::array mCache{{}}; ///< caching positions, covariances and angles for track extrapolations and interpolation std::vector mGIDsSuccess; ///< keep track of the GIDs which could be processed successfully + TrackValidationData mTrackValidation; + // helpers o2::gpu::GPUTRDRecoParam mRecoParam; ///< parameters required for TRD refit o2::trd::Geometry* mGeoTRD; ///< TRD geometry instance (needed for tilted pad correction) @@ -447,6 +489,9 @@ class TrackInterpolation float mBz; ///< required for helix approximation bool mInitDone{false}; ///< initialization done flag size_t mRejectedResiduals{}; ///< number of rejected residuals + size_t mNRejRefit = 0; + size_t mNRejProp = 0; + size_t mNRejLoop = 0; ClassDefNV(TrackInterpolation, 1); }; diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h index c9226589ec703..202df063bf0fc 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h @@ -259,12 +259,13 @@ class TrackResiduals /// \param r fit result for circle radius is stored here /// \param residHelixY residuals in y from fitted circle to given points is stored here static void fitCircle(int nCl, std::array& x, std::array& y, float& xc, float& yc, float& r, std::array& residHelixY); - + static void fitCircle(TrackInterpolation::TrackValidationData& params); /// Fits a straight line to a given set of points, w/o taking into account measurement errors or different weights for the points /// Straight line is given by y = a * x + b /// \param res[0] contains the slope (a) /// \param res[1] contains the offset (b) static bool fitPoly1(int nCl, std::array& x, std::array& y, std::array& res); + static bool fitPoly1(TrackInterpolation::TrackValidationData& params); // -------------------------------------- binning / geometry -------------------------------------------------- diff --git a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h index 4703c7ff39fce..e77610acb8e7e 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h +++ b/Detectors/TPC/calibration/SpacePoints/src/SpacePointCalibLinkDef.h @@ -38,6 +38,9 @@ #pragma link C++ class o2::calibration::TimeSlot < o2::tpc::ResidualsContainer> + ; #pragma link C++ class o2::calibration::TimeSlotCalibration < o2::tpc::ResidualsContainer> + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::tpc::SpacePointsCalibConfParam> + ; +#pragma link C++ class o2::tpc::TrackInterpolation::ValidationPoint + ; +#pragma link C++ class std::vector < o2::tpc::TrackInterpolation::ValidationPoint> + ; +#pragma link C++ class o2::tpc::TrackInterpolation::TrackValidationData + ; #pragma link C++ struct o2::tpc::SpacePointsCalibConfParam; #pragma read sourceClass = "o2::tpc::TrackData" targetClass = "o2::tpc::TrackData" source = "o2::track::TrackPar par" version = "[-10]" target = "par" code = "{}"; diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index 76daab93dd8e0..abb06292fe25d 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -111,6 +111,19 @@ void UnbinnedResid::init(long timestamp) gInitDone = true; } +TrackInterpolation::~TrackInterpolation() +{ + finalize(); +} + +void TrackInterpolation::finalize() +{ + if (mDBGOut) { + mDBGOut->Close(); + mDBGOut.reset(); + } +} + void TrackInterpolation::init(o2::dataformats::GlobalTrackID::mask_t src, o2::dataformats::GlobalTrackID::mask_t srcMap) { // perform initialization @@ -141,6 +154,12 @@ void TrackInterpolation::init(o2::dataformats::GlobalTrackID::mask_t src, o2::da auto geom = o2::its::GeometryTGeo::Instance(); geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); mTPCParam = o2::gpu::GPUO2InterfaceUtils::getFullParamShared(0.f, mNHBPerTF); + + if (mParams->writeValidationData) { + std::string dbgnm = mNLanes == 1 ? "track_interpolation_dbg.root" : fmt::format("track_interpolation_dbg_{}.root", mLaneID); + mDBGOut = std::make_unique(dbgnm.c_str(), "recreate"); + } + mInitDone = true; LOGP(info, "Done initializing TrackInterpolation. Configured track input: {}. Track input specifically for map: {}", GTrackID::getSourcesNames(mSourcesConfigured), mSingleSourcesConfigured ? "identical" : GTrackID::getSourcesNames(mSourcesConfiguredMap)); @@ -252,7 +271,7 @@ void TrackInterpolation::prepareInputTrackSample(const o2::globaltracking::RecoC } } - for (int is = GTrackID::NSources; is >= 0; is--) { + for (int is = GTrackID::NSources; is--;) { if (!allowedSources[is]) { continue; } @@ -355,7 +374,6 @@ void TrackInterpolation::process() trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTRD]].end()); trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPCTOF]].end()); trackIndices.insert(trackIndices.end(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].begin(), mTrackIndices[mTrackTypes[GTrackID::ITSTPC]].end()); - int nSeeds = mSeeds.size(), lastChecked = 0; mParentID.clear(); mParentID.resize(nSeeds, -1); @@ -420,7 +438,9 @@ void TrackInterpolation::process() remSeeds.resize(mSeeds.size() - lastChecked); std::iota(remSeeds.begin(), remSeeds.end(), lastChecked); std::shuffle(remSeeds.begin(), remSeeds.end(), g); - LOGP(info, "Up to {} tracks out of {} additional seeds will be processed in random order, of which {} are stripped versions, accepted seeds: {}", mAddTracksForMapPerTF, remSeeds.size(), mSeeds.size() - nSeeds, mTrackDataCompact.size()); + LOGP(info, "Up to {} tracks out of {} additional seeds will be processed in random order, of which {} are stripped versions, accepted seeds: {}", + mAddTracksForMapPerTF > 0 ? mAddTracksForMapPerTF : remSeeds.size(), + remSeeds.size(), mSeeds.size() - nSeeds, mTrackDataCompact.size()); } int extraChecked = 0; for (int iSeed : remSeeds) { @@ -437,8 +457,12 @@ void TrackInterpolation::process() extrapolateTrack(iSeed); } } - LOG(info) << "Could process " << mTrackData.size() << " tracks successfully. " << mRejectedResiduals << " residuals were rejected. " << mClRes.size() << " residuals were accepted."; + LOGP(info, "Could process {} tracks successfully ({} rejected in refits, {} in propagation, {} as loopers), {} residuals were rejected, {} accepted", + mTrackData.size(), mNRejRefit, mNRejProp, mNRejLoop, mRejectedResiduals, mClRes.size()); mRejectedResiduals = 0; + mNRejRefit = 0; + mNRejProp = 0; + mNRejLoop = 0; } void TrackInterpolation::interpolateTrack(int iSeed) @@ -467,6 +491,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } } if (mParams->refitITS && !refITSTrack(gidTable[GTrackID::ITS], iSeed)) { + mNRejRefit++; return; } trackData.gid = mGIDs[iSeed]; @@ -510,10 +535,12 @@ void TrackInterpolation::interpolateTrack(int iSeed) } if (!trkWork.rotate(mCache[iRow].clAngle)) { LOG(debug) << "Failed to rotate track during first extrapolation"; + mNRejProp++; return; } if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on first extrapolation"; + mNRejProp++; return; } mCache[iRow].y[ExtOut] = trkWork.getY(); @@ -537,6 +564,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) const float clTOFAlpha = o2::math_utils::sector2Angle(clTOFSec); if (!trkWork.rotate(clTOFAlpha)) { LOG(debug) << "Failed to rotate into TOF cluster sector frame"; + mNRejProp++; return; } float clTOFxyz[3] = {clTOF.getX(), clTOF.getY(), clTOF.getZ()}; @@ -547,12 +575,14 @@ void TrackInterpolation::interpolateTrack(int iSeed) std::array clTOFCov{mParams->sigYZ2TOF, 0.f, mParams->sigYZ2TOF}; // assume no correlation between y and z and equal cluster error sigma^2 = (3cm)^2 / 12 if (!propagator->PropagateToXBxByBz(trkWork, clTOFxyz[0], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed final propagation to TOF radius"; + mNRejProp++; return; } // TODO: check if reset of covariance matrix is needed here (or, in case TOF point is not available at outermost TRD layer) if (!trkWork.update(clTOFYZ, clTOFCov)) { LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster"; // LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]); + mNRejProp++; return; } } @@ -574,6 +604,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) } if (!trkWork.update(trkltTRDYZ, trkltTRDCov)) { LOG(debug) << "Failed to update track at TRD layer " << iLayer; + mNRejProp++; return; } } @@ -601,11 +632,13 @@ void TrackInterpolation::interpolateTrack(int iSeed) } if (!trkWork.rotate(mCache[iRow].clAngle)) { LOG(debug) << "Failed to rotate track during back propagation"; + mNRejProp++; return; } if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) { LOG(debug) << "Failed on back propagation"; // printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha); + mNRejProp++; return; } mCache[iRow].y[ExtIn] = trkWork.getY(); @@ -665,15 +698,17 @@ void TrackInterpolation::interpolateTrack(int iSeed) } trackData.dEdxTPC = trkTPC.getdEdx().dEdxTotTPC; - TrackParams params; // for refitted track parameters and flagging rejected clusters - if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) { - // track is good + mTrackValidation.clear(); // for refitted track parameters and flagging rejected clusters + + bool stored = false; + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals); + if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0; int iRow = 0; for (unsigned int iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; - if (params.flagRej[iCl]) { - // skip masked cluster residual + const auto rej = trackData.filterFlag < 0 ? false : mTrackValidation.points[iCl].flagRej; + if (rej && !mParams->keepRejectedResiduals) { // skip masked cluster residual continue; } const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp)); @@ -683,7 +718,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) const auto z = clusterResiduals[iCl].z; const auto sec = clusterResiduals[iCl].sec; if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(y) < param::MaxY) && (std::abs(z) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) { - mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, sec); + mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, sec, -1, rej); mDetInfoRes.emplace_back().setTPC(mCacheDEDX[iRow].first, mCacheDEDX[iRow].second); // qtot, qmax ++nClValidated; } else { @@ -836,20 +871,18 @@ void TrackInterpolation::interpolateTrack(int iSeed) } mGIDsSuccess.push_back(mGIDs[iSeed]); - mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid); + mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); mTrackData.push_back(std::move(trackData)); + stored = true; if (mDumpTrackPoints) { (*trackDataExtended).clIdx.setEntries(nClValidated); (*trackDataExtended).nExtDetResid = trackData.nExtDetResid; + (*trackDataExtended).filterFlag = trackData.filterFlag; mTrackDataExtended.push_back(std::move(*trackDataExtended)); } } - if (mParams->writeUnfiltered) { - TrackData trkDataTmp = trackData; - trkDataTmp.clIdx.setFirstEntry(mClResUnfiltered.size()); - trkDataTmp.clIdx.setEntries(clusterResiduals.size()); - mTrackDataUnfiltered.push_back(std::move(trkDataTmp)); - mClResUnfiltered.insert(mClResUnfiltered.end(), clusterResiduals.begin(), clusterResiduals.end()); + if (mParams->writeValidationData && trackData.filterFlag >= 0 && mDBGOut) { + (*mDBGOut) << "valdata" << "params=" << mTrackValidation << "trackData=" << (stored ? mTrackData.back() : trackData) << "\n"; } } @@ -882,7 +915,7 @@ int TrackInterpolation::processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLa float tiltCorrUp = tilt * (trdSP.getZ() - trkWork.getZ()); float zPosCorrUp = trdSP.getZ() + mRecoParam.getZCorrCoeffNRC() * trkWork.getTgl(); // maybe Z can be corrected on avarage already by the tracklet transformer? float padLength = pad->getRowSize(trdTrklt.getPadRow()); - if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(trdSP.getZ() - trkWork.getZ()) < padLength))) { + if (!((trkWork.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::abs(trdSP.getZ() - trkWork.getZ()) < padLength))) { tiltCorrUp = 0.f; } (*trkltTRDYZ)[0] = trdSP.getY() - tiltCorrUp; @@ -933,6 +966,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } } if (mParams->refitITS && !refITSTrack(gidTable[GTrackID::ITS], iSeed)) { + mNRejRefit++; return; } trackData.gid = mGIDs[iSeed]; @@ -957,6 +991,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) // we seem to be looping, abort this track LOGP(debug, "TPC track with pT={} GeV and {} clusters has cluster {} on row {} while the previous cluster was on row {}", mSeeds[iSeed].getPt(), trkTPC.getNClusterReferences(), iCl, row, clRowPrev); + mNRejLoop++; return; } else { // this is the first cluster we see on this pad row @@ -965,9 +1000,11 @@ void TrackInterpolation::extrapolateTrack(int iSeed) float x = 0, y = 0, z = 0; mFastTransform->TransformIdeal(sector, row, cl.getPad(), cl.getTime(), x, y, z, clusterTimeBinOffset); if (!trkWork.rotate(o2::math_utils::sector2Angle(sector))) { + mNRejProp++; return; } if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) { + mNRejProp++; return; } @@ -988,9 +1025,10 @@ void TrackInterpolation::extrapolateTrack(int iSeed) ++nMeasurements; } - TrackParams params; // for refitted track parameters and flagging rejected clusters + mTrackValidation.clear(); // for refitted track parameters and flagging rejected clusters if (clusterResiduals.size() > constants::MAXGLOBALPADROW) { LOGP(warn, "Extrapolated ITS-TPC track and found more residuals than possible ({})", clusterResiduals.size()); + mNRejLoop++; return; } @@ -1004,14 +1042,18 @@ void TrackInterpolation::extrapolateTrack(int iSeed) (*trackDataExtended).trkOuter = trkWork; } - if (mParams->skipOutlierFiltering || validateTrack(trackData, params, clusterResiduals)) { - // track is good, store TPC part - + bool stored = false; + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals); + if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0, iRow = 0; unsigned int iCl = 0; for (iCl = 0; iCl < clusterResiduals.size(); ++iCl) { iRow += clusterResiduals[iCl].dRow; - if (iRow < param::NPadRows && params.flagRej[iCl]) { // skip masked cluster residual + if (iRow >= param::NPadRows) { // RS why do we need this? + continue; + } + const auto rej = trackData.filterFlag < 0 ? false : mTrackValidation.points[iCl].flagRej; + if (rej && !mParams->keepRejectedResiduals) { // skip masked cluster residual continue; } const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp)); @@ -1020,7 +1062,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) const auto y = clusterResiduals[iCl].y; const auto z = clusterResiduals[iCl].z; if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(y) < param::MaxY) && (std::abs(z) < param::MaxZ) && (std::abs(tgPhi) < param::MaxTgSlp)) { - mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, clusterResiduals[iCl].sec); + mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, clusterResiduals[iCl].sec, -1, rej); mDetInfoRes.emplace_back().setTPC(mCacheDEDX[iRow].first, mCacheDEDX[iRow].second); // qtot, qmax ++nClValidated; } else { @@ -1183,88 +1225,95 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } } mTrackData.push_back(std::move(trackData)); + stored = true; mGIDsSuccess.push_back(mGIDs[iSeed]); - mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid); + mTrackDataCompact.emplace_back(trackData.clIdx.getFirstEntry(), trackData.multStack, nClValidated, mGIDs[iSeed].getSource(), trackData.nExtDetResid, trackData.filterFlag); if (mDumpTrackPoints) { (*trackDataExtended).clIdx.setEntries(nClValidated); (*trackDataExtended).nExtDetResid = trackData.nExtDetResid; + (*trackDataExtended).filterFlag = trackData.filterFlag; mTrackDataExtended.push_back(std::move(*trackDataExtended)); } } - if (mParams->writeUnfiltered) { - TrackData trkDataTmp = trackData; - trkDataTmp.clIdx.setFirstEntry(mClResUnfiltered.size()); - trkDataTmp.clIdx.setEntries(clusterResiduals.size()); - mTrackDataUnfiltered.push_back(std::move(trkDataTmp)); - mClResUnfiltered.insert(mClResUnfiltered.end(), clusterResiduals.begin(), clusterResiduals.end()); + if (mParams->writeValidationData && trackData.filterFlag >= 0 && mDBGOut) { + (*mDBGOut) << "valdata" << "params=" << mTrackValidation << "trackData=" << (stored ? mTrackData.back() : trackData) << "\n"; } } -bool TrackInterpolation::validateTrack(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const +int8_t TrackInterpolation::validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes) { - if (clsRes.size() < mParams->minNCl) { - // no enough clusters for this track to be considered - LOG(debug) << "Skipping track with too few clusters: " << clsRes.size(); - return false; - } + int8_t status = 0; + while (true) { + if (clsRes.size() < mParams->minNCl) { + // no enough clusters for this track to be considered + LOG(debug) << "Skipping track with too few clusters: " << clsRes.size(); + status |= 0x1; + if (!mParams->keepRejectedResiduals) { + break; // we don't keep de-validated tracks, no need to check further + } + } - bool resHelix = compareToHelix(trk, params, clsRes); - if (!resHelix) { - LOG(debug) << "Skipping track too far from helix approximation"; - return false; - } - if (fabsf(mBz) > 0.01 && fabsf(params.qpt) > mParams->maxQ2Pt) { - LOG(debug) << "Skipping track with too high q/pT: " << params.qpt; - return false; - } - if (!outlierFiltering(trk, params, clsRes)) { - return false; + bool resHelix = compareToHelix(trk, params, clsRes); + if (!resHelix) { + LOG(debug) << "Skipping track too far from helix approximation"; + status |= 0x1 << 1; + if (!mParams->keepRejectedResiduals) { + break; // we don't keep de-validated tracks, no need to check further + } + } + if (std::abs(mBz) > 0.01 && std::abs(params.qpt) > mParams->maxQ2Pt) { + LOG(debug) << "Skipping track with too high q/pT: " << params.qpt; + status |= 0x1 << 2; + if (!mParams->keepRejectedResiduals) { + break; // we don't keep de-validated tracks, no need to check further + } + } + if (!outlierFiltering(trk, params, clsRes)) { + status |= 0x1 << 3; + if (!mParams->keepRejectedResiduals) { + break; // we don't keep de-validated tracks, no need to check further + } + } + break; } - return true; + return status & 0x7f; } -bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const +bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes) { - std::array residHelixY; - std::array residHelixZ; - - std::array xLab; - std::array yLab; - std::array sPath; - - float curvature = fabsf(trk.par.getQ2Pt() * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f); + float curvature = std::abs(trk.par.getQ2Pt() * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f); int secFirst = clsRes[0].sec; float phiSect = (secFirst + .5f) * o2::constants::math::SectorSpanRad; float snPhi = sin(phiSect); float csPhi = cos(phiSect); - sPath[0] = 0.f; int iRow = 0; int nCl = clsRes.size(); for (unsigned int iP = 0; iP < nCl; ++iP) { + auto& point = params.points.emplace_back(); + iRow += clsRes[iP].dRow; - float yTrk = clsRes[iP].y; - // LOGF(info, "iRow(%i), yTrk(%f)", iRow, yTrk); - xLab[iP] = param::RowX[iRow]; + point.yTrk = clsRes[iP].y; + point.sec = clsRes[iP].sec; if (clsRes[iP].sec != secFirst) { float phiSectCurrent = (clsRes[iP].sec + .5f) * o2::constants::math::SectorSpanRad; float cs = cos(phiSectCurrent - phiSect); float sn = sin(phiSectCurrent - phiSect); - xLab[iP] = param::RowX[iRow] * cs - yTrk * sn; - yLab[iP] = yTrk * cs + param::RowX[iRow] * sn; + point.xLab = param::RowX[iRow] * cs - point.yTrk * sn; + point.yLab = point.yTrk * cs + param::RowX[iRow] * sn; } else { - xLab[iP] = param::RowX[iRow]; - yLab[iP] = yTrk; + point.xLab = param::RowX[iRow]; + point.yLab = point.yTrk; } // this is needed only later, but we retrieve it already now to save another loop - params.zTrk[iP] = clsRes[iP].z; - params.xTrk[iP] = param::RowX[iRow]; - params.dy[iP] = clsRes[iP].dy; - params.dz[iP] = clsRes[iP].dz; + point.zTrk = clsRes[iP].z; + point.xTrk = param::RowX[iRow]; + point.dy = clsRes[iP].dy; + point.dz = clsRes[iP].dz; // done retrieving values for later if (iP > 0) { - float dx = xLab[iP] - xLab[iP - 1]; - float dy = yLab[iP] - yLab[iP - 1]; + float dx = point.xLab - params.points[iP - 1].xLab; + float dy = point.yLab - params.points[iP - 1].yLab; float ds2 = dx * dx + dy * dy; float ds = sqrt(ds2); // circular path (linear approximation) // if the curvature of the track or the (approximated) chord length is too large the more exact formula is used: @@ -1273,26 +1322,19 @@ bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& param if (ds * curvature > 0.05) { ds *= (1.f + ds2 * curvature * curvature / 24.f); } - sPath[iP] = sPath[iP - 1] + ds; + point.sPath = params.points[iP - 1].sPath + ds; + } else { + point.sPath = 0; } } - if (fabsf(mBz) < 0.01) { + if (std::abs(mBz) < 0.01) { // for B=0 we don't need to try a circular fit... return true; } - float xcSec = 0.f; - float ycSec = 0.f; - float r = 0.f; - TrackResiduals::fitCircle(nCl, xLab, yLab, xcSec, ycSec, r, residHelixY); - // LOGF(info, "Done with circle fit. nCl(%i), xcSec(%f), ycSec(%f), r(%f).", nCl, xcSec, ycSec, r); - /* - for (int i=0; i pol1Z; - TrackResiduals::fitPoly1(nCl, sPath, params.zTrk, pol1Z); + params.qpt = std::copysign(1.f / (params.r * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f), curvSign); - params.tgl = pol1Z[0]; + TrackResiduals::fitPoly1(params); // max deviations in both directions from helix fit in y and z float hMinY = 1e9f; @@ -1325,23 +1360,24 @@ bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& param float hMinZ = 1e9f; float hMaxZ = -1e9f; // extract residuals in Z and fill track slopes in sector frame - int secCurr = secFirst; + int secCurr = -1; iRow = 0; + float xcSec = 0; for (unsigned int iCl = 0; iCl < nCl; ++iCl) { iRow += clsRes[iCl].dRow; - float resZ = params.zTrk[iCl] - (pol1Z[1] + sPath[iCl] * pol1Z[0]); - residHelixZ[iCl] = resZ; - if (resZ < hMinZ) { - hMinZ = resZ; + auto& pnt = params.points[iCl]; + pnt.residHelixZ = pnt.zTrk - (params.zOffs + pnt.sPath * params.tgl); + if (pnt.residHelixZ < hMinZ) { + hMinZ = pnt.residHelixZ; } - if (resZ > hMaxZ) { - hMaxZ = resZ; + if (pnt.residHelixZ > hMaxZ) { + hMaxZ = pnt.residHelixZ; } - if (residHelixY[iCl] < hMinY) { - hMinY = residHelixY[iCl]; + if (pnt.residHelixY < hMinY) { + hMinY = pnt.residHelixY; } - if (residHelixY[iCl] > hMaxY) { - hMaxY = residHelixY[iCl]; + if (pnt.residHelixY > hMaxY) { + hMaxY = pnt.residHelixY; } int sec = clsRes[iCl].sec; if (sec != secCurr) { @@ -1349,34 +1385,34 @@ bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& param phiSect = (.5f + sec) * o2::constants::math::SectorSpanRad; snPhi = sin(phiSect); csPhi = cos(phiSect); - xcSec = xc * csPhi + yc * snPhi; // recalculate circle center in the sector frame + xcSec = params.xcLab * csPhi + params.ycLab * snPhi; // recalculate circle center in the new sector frame } - float cstalp = (param::RowX[iRow] - xcSec) / r; - if (fabsf(cstalp) > 1.f - sFloatEps) { + float cstalp = (param::RowX[iRow] - xcSec) / params.r; + if (std::abs(cstalp) > 1.f - sFloatEps) { // track cannot reach this pad row cstalp = std::copysign(1.f - sFloatEps, cstalp); } - params.tglArr[iCl] = cstalp / sqrt((1 - cstalp) * (1 + cstalp)); // 1 / tan(acos(cstalp)) = cstalp / sqrt(1 - cstalp^2) + pnt.tglArr = cstalp / sqrt((1 - cstalp) * (1 + cstalp)); // 1 / tan(acos(cstalp)) = cstalp / sqrt(1 - cstalp^2) // In B+ the slope of q- should increase with x. Just look on q * B if (params.qpt * mBz > 0) { - params.tglArr[iCl] *= -1.f; + pnt.tglArr = -pnt.tglArr; } } // LOGF(info, "CompareToHelix: hMaxY(%f), hMinY(%f), hMaxZ(%f), hMinZ(%f). Max deviation allowed: y(%.2f), z(%.2f)", hMaxY, hMinY, hMaxZ, hMinZ, mParams->maxDevHelixY, mParams->maxDevHelixZ); // LOGF(info, "New pt/Q (%f), old pt/Q (%f)", 1./params.qpt, 1./trk.qPt); - return fabsf(hMaxY - hMinY) < mParams->maxDevHelixY && fabsf(hMaxZ - hMinZ) < mParams->maxDevHelixZ; + return std::abs(hMaxY - hMinY) < mParams->maxDevHelixY && std::abs(hMaxZ - hMinZ) < mParams->maxDevHelixZ; } -bool TrackInterpolation::outlierFiltering(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const +bool TrackInterpolation::outlierFiltering(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes) { if (clsRes.size() < mParams->nMALong) { LOG(debug) << "Skipping track with too few clusters for long moving average: " << clsRes.size(); return false; } float rmsLong = checkResiduals(trk, params, clsRes); - if (static_cast(params.flagRej.count()) / clsRes.size() > mParams->maxRejFrac) { - LOGP(debug, "Skipping track with too many clusters rejected: {} out of {}", params.flagRej.count(), clsRes.size()); + if (static_cast(params.nRej) / clsRes.size() > mParams->maxRejFrac) { + LOGP(debug, "Skipping track with too many clusters rejected: {} out of {}", params.nRej, clsRes.size()); return false; } if (rmsLong > mParams->maxRMSLong) { @@ -1386,7 +1422,7 @@ bool TrackInterpolation::outlierFiltering(const TrackData& trk, TrackParams& par return true; } -float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& params, const std::vector& clsRes) const +float TrackInterpolation::checkResiduals(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes) { float rmsLong = 0.f; @@ -1395,9 +1431,14 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para int iClLast = nCl - 1; int secStart = clsRes[0].sec; + auto rejectAll = [¶ms]() { + for (auto& pnt : params.points) { + pnt.flagRej = true; + } + params.nRej = params.points.size(); + }; + // arrays with differences / abs(differences) of points to their neighbourhood, initialized to zero - std::array yDiffLL{}; - std::array zDiffLL{}; std::array absDevY{}; std::array absDevZ{}; @@ -1411,8 +1452,7 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para if (iCl == iClLast) { ++nClSec; } - diffToLocLine(nClSec, iClFirst, params.xTrk, params.dy, yDiffLL); - diffToLocLine(nClSec, iClFirst, params.xTrk, params.dz, zDiffLL); + diffToLocLine(params, iClFirst, nClSec); iClFirst = iCl; secStart = clsRes[iCl].sec; } @@ -1420,17 +1460,18 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para int nAccY = 0; int nAccZ = 0; for (int iCl = nCl; iCl--;) { - if (fabsf(yDiffLL[iCl]) > param::sEps) { - absDevY[nAccY++] = fabsf(yDiffLL[iCl]); + const auto pnt = params.points[iCl]; + if (std::abs(pnt.diffYSmooth) > param::sEps) { + absDevY[nAccY++] = std::abs(pnt.diffYSmooth); } - if (fabsf(zDiffLL[iCl]) > param::sEps) { - absDevZ[nAccZ++] = fabsf(zDiffLL[iCl]); + if (std::abs(pnt.diffZSmooth) > param::sEps) { + absDevZ[nAccZ++] = std::abs(pnt.diffZSmooth); } } if (nAccY < mParams->minNumberOfAcceptedResiduals || nAccZ < mParams->minNumberOfAcceptedResiduals) { // mask all clusters LOGP(debug, "Accepted {} clusters for dY {} clusters for dZ, but required at least {} for both", nAccY, nAccZ, mParams->minNumberOfAcceptedResiduals); - params.flagRej.set(); + rejectAll(); return 0.f; } // estimate rms on 90% of the smallest deviations @@ -1450,7 +1491,7 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para rmsZkeep = std::sqrt(rmsZkeep / nKeepZ); if (rmsYkeep < param::sEps || rmsZkeep < param::sEps) { LOG(warning) << "Too small RMS: " << rmsYkeep << "(y), " << rmsZkeep << "(z)."; - params.flagRej.set(); + rejectAll(); return 0.f; } float rmsYkeepI = 1.f / rmsYkeep; @@ -1459,18 +1500,19 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para std::array yAcc; std::array yDiffLong; for (int iCl = 0; iCl < nCl; ++iCl) { - yDiffLL[iCl] *= rmsYkeepI; - zDiffLL[iCl] *= rmsZkeepI; - if (yDiffLL[iCl] * yDiffLL[iCl] + zDiffLL[iCl] * zDiffLL[iCl] > mParams->maxStdDevMA) { - params.flagRej.set(iCl); + auto& pnt = params.points[iCl]; + auto yDiffScl = pnt.diffYSmooth * rmsYkeepI; + auto zDiffScl = pnt.diffZSmooth * rmsZkeepI; + if (yDiffScl * yDiffScl + zDiffScl * zDiffScl > mParams->maxStdDevMA) { + pnt.flagRej = true; + params.nRej++; } else { - yAcc[nAcc++] = params.dy[iCl]; + yAcc[nAcc++] = pnt.dy; } } if (nAcc > mParams->nMALong) { diffToMA(nAcc, yAcc, yDiffLong); - float average = 0.f; - float rms = 0.f; + float average = 0.f, rms = 0.f; for (int i = 0; i < nAcc; ++i) { average += yDiffLong[i]; rms += yDiffLong[i] * yDiffLong[i]; @@ -1482,86 +1524,72 @@ float TrackInterpolation::checkResiduals(const TrackData& trk, TrackParams& para return rmsLong; } -void TrackInterpolation::diffToLocLine(const int np, int idxOffset, const std::array& x, const std::array& y, std::array& diffY) const +void TrackInterpolation::diffToLocLine(TrackValidationData& params, int start, int np) { - // Calculate the difference between the points and the linear extrapolations from the neighbourhood. - // Nothing more than multiple 1-d fits at once. Instead of building 4 sums (x, x^2, y, xy), 4 * nPoints sums are calculated at once - // compare to TrackResiduals::fitPoly1() method - - // adding one entry to the vectors saves an additional if statement when calculating the cumulants - std::vector sumX1vec(np + 1); - std::vector sumX2vec(np + 1); - std::vector sumY1vec(np + 1); - std::vector sumXYvec(np + 1); - auto sumX1 = &(sumX1vec[1]); - auto sumX2 = &(sumX2vec[1]); - auto sumY1 = &(sumY1vec[1]); - auto sumXY = &(sumXYvec[1]); - - // accumulate sums for all points - for (int iCl = 0; iCl < np; ++iCl) { - int idx = iCl + idxOffset; - sumX1[iCl] = sumX1[iCl - 1] + x[idx]; - sumX2[iCl] = sumX2[iCl - 1] + x[idx] * x[idx]; - sumY1[iCl] = sumY1[iCl - 1] + y[idx]; - sumXY[iCl] = sumXY[iCl - 1] + x[idx] * y[idx]; - } - - for (int iCl = 0; iCl < np; ++iCl) { - int iClLeft = iCl - mParams->nMAShort; - int iClRight = iCl + mParams->nMAShort; - if (iClLeft < 0) { - iClLeft = 0; - } - if (iClRight >= np) { - iClRight = np - 1; - } - int nPoints = iClRight - iClLeft; + std::array sumX1{}, sumX2{}, sumY1{}, sumXY{}, sumZ1{}, sumXZ{}; + for (int i = 0; i < np; ++i) { + const auto& pnt = params.points[start + i]; + const float x = pnt.xTrk, y = pnt.dy, z = pnt.dz; + sumX1[i + 1] = sumX1[i] + x; + sumX2[i + 1] = sumX2[i] + x * x; + sumY1[i + 1] = sumY1[i] + y; + sumXY[i + 1] = sumXY[i] + x * y; + sumZ1[i + 1] = sumZ1[i] + z; + sumXZ[i + 1] = sumXZ[i] + x * z; + } + + for (int i = 0; i < np; ++i) { + auto& pnt = params.points[start + i]; + + const int iLeft = std::max(0, i - mParams->nMAShort); + const int iRight = std::min(np - 1, i + mParams->nMAShort); + + const int nPoints = iRight - iLeft; // excluding current point + if (nPoints < mParams->nMAShort) { continue; } - float nPointsInv = 1.f / nPoints; - int iClLeftP = iClLeft - 1; - int iClCurrP = iCl - 1; - // extract sum from iClLeft to iClRight from cumulants, excluding iCl from the fit - float sX1 = sumX1[iClRight] - sumX1[iClLeftP] - (sumX1[iCl] - sumX1[iClCurrP]); - float sX2 = sumX2[iClRight] - sumX2[iClLeftP] - (sumX2[iCl] - sumX2[iClCurrP]); - float sY1 = sumY1[iClRight] - sumY1[iClLeftP] - (sumY1[iCl] - sumY1[iClCurrP]); - float sXY = sumXY[iClRight] - sumXY[iClLeftP] - (sumXY[iCl] - sumXY[iClCurrP]); - float det = sX2 - nPointsInv * sX1 * sX1; - if (fabsf(det) < 1e-12f) { + + const float nPointsInv = 1.f / nPoints; + + float sX1 = sumX1[iRight + 1] - sumX1[iLeft] - pnt.xTrk; + float sX2 = sumX2[iRight + 1] - sumX2[iLeft] - pnt.xTrk * pnt.xTrk; + float sY1 = sumY1[iRight + 1] - sumY1[iLeft] - pnt.dy; + float sXY = sumXY[iRight + 1] - sumXY[iLeft] - pnt.xTrk * pnt.dy; + float sZ1 = sumZ1[iRight + 1] - sumZ1[iLeft] - pnt.dz; + float sXZ = sumXZ[iRight + 1] - sumXZ[iLeft] - pnt.xTrk * pnt.dz; + + const float det = sX2 - nPointsInv * sX1 * sX1; + + if (std::abs(det) < 1e-12f) { continue; } - float slope = (sXY - nPointsInv * sX1 * sY1) / det; - float offset = nPointsInv * sY1 - nPointsInv * slope * sX1; - diffY[iCl + idxOffset] = y[iCl + idxOffset] - slope * x[iCl + idxOffset] - offset; + + const float slopeY = (sXY - nPointsInv * sX1 * sY1) / det; + const float offsetY = nPointsInv * (sY1 - slopeY * sX1); + const float slopeZ = (sXZ - nPointsInv * sX1 * sZ1) / det; + const float offsetZ = nPointsInv * (sZ1 - slopeZ * sX1); + pnt.diffYSmooth = pnt.dy - (slopeY * pnt.xTrk + offsetY); + pnt.diffZSmooth = pnt.dz - (slopeZ * pnt.xTrk + offsetZ); } } -void TrackInterpolation::diffToMA(const int np, const std::array& y, std::array& diffMA) const +void TrackInterpolation::diffToMA(const int np, const std::array& y, std::array& diffMA) { // Calculate - std::vector sumVec(np + 1); - auto sum = &(sumVec[1]); + std::array sum{}; for (int i = 0; i < np; ++i) { - sum[i] = sum[i - 1] + y[i]; + sum[i + 1] = sum[i] + y[i]; } for (int i = 0; i < np; ++i) { - diffMA[i] = 0; - int iLeft = i - mParams->nMALong; - int iRight = i + mParams->nMALong; - if (iLeft < 0) { - iLeft = 0; - } - if (iRight >= np) { - iRight = np - 1; - } + diffMA[i] = 0.f; + int iLeft = std::max(0, i - mParams->nMALong); + int iRight = std::min(np - 1, i + mParams->nMALong); int nPoints = iRight - iLeft; - if (nPoints < mParams->nMALong) { - // this cannot happen, since at least mParams->nMALong points are required as neighbours for this function to be called + if (nPoints < mParams->nMALong) { // this cannot happen, since at least mParams->nMALong points are required as neighbours for this function to be called continue; } - float movingAverage = (sum[iRight] - sum[iLeft - 1] - (sum[i] - sum[i - 1])) / nPoints; + float movingAverage = (sum[iRight + 1] - sum[iLeft] - y[i]) / nPoints; diffMA[i] = y[i] - movingAverage; } } @@ -1573,8 +1601,6 @@ void TrackInterpolation::reset() mTrackDataExtended.clear(); mClRes.clear(); mDetInfoRes.clear(); - mTrackDataUnfiltered.clear(); - mClResUnfiltered.clear(); mGIDsSuccess.clear(); for (auto& vec : mTrackIndices) { vec.clear(); diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx index d3db11daf9e87..49528f581d2b6 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx @@ -84,7 +84,7 @@ void TrackResiduals::setY2XBinning(const std::vector& binning) } const int nBins = binning.size() - 1; - if (fabsf(binning[0] + 1.f) > param::sEps || fabsf(binning[nBins] - 1.f) > param::sEps) { + if (std::abs(binning[0] + 1.f) > param::sEps || std::abs(binning[nBins] - 1.f) > param::sEps) { LOG(error) << "Provided binning for y/x not in range -1 to 1: " << binning[0] << " - " << binning[nBins] << ". Not changing y/x binning"; return; } @@ -122,7 +122,7 @@ void TrackResiduals::setZ2XBinning(const std::vector& binning) } int nBins = binning.size() - 1; - if (fabsf(binning[0]) > param::sEps || fabsf(binning[nBins] - 1.f) > param::sEps) { + if (std::abs(binning[0]) > param::sEps || std::abs(binning[nBins] - 1.f) > param::sEps) { LOG(error) << "Provided binning for z/x not in range 0 to 1: " << binning[0] << " - " << binning[nBins] << ". Not changing z/x binning"; return; } @@ -267,7 +267,7 @@ int TrackResiduals::getRowID(float x) const bool TrackResiduals::findVoxelBin(int secID, float x, float y, float z, std::array& bvox) const { // Z/X bin - if (fabs(z / x) > mMaxZ2X) { + if (std::abs(z / x) > mMaxZ2X) { return false; } int bz = getZ2XBinExact(secID < SECTORSPERSIDE ? z / x : -z / x); @@ -585,7 +585,7 @@ int TrackResiduals::validateVoxels(int iSec) // check fit errors if (resVox.E[ResY] * resVox.E[ResY] > mParams->maxFitErrY2 || resVox.E[ResX] * resVox.E[ResX] > mParams->maxFitErrX2 || - fabs(resVox.EXYCorr) > mParams->maxFitCorrXY) { + std::abs(resVox.EXYCorr) > mParams->maxFitCorrXY) { voxelOK = false; ++cntMaskedFit; } @@ -973,7 +973,7 @@ bool TrackResiduals::getSmoothEstimate(int iSec, float x, float p, float z, std: } double vi = voxNb->D[iDim]; double wi = wiCache; - if (mUseErrInSmoothing && fabs(voxNb->E[iDim]) > 1e-6) { + if (mUseErrInSmoothing && std::abs(voxNb->E[iDim]) > 1e-6) { // account for point error apart from kernel value wi /= (voxNb->E[iDim] * voxNb->E[iDim]); } @@ -1222,7 +1222,7 @@ void TrackResiduals::medFit(int nPoints, int offset, const std::vector& x if (sigb > 0) { float b2 = bb + std::copysign(3.f * sigb, f1); float f2 = roFunc(nPoints, offset, x, y, b2, aa); - if (fabs(f1 - f2) < sFloatEps) { + if (std::abs(f1 - f2) < sFloatEps) { a = aa; b = bb; return; @@ -1235,7 +1235,7 @@ void TrackResiduals::medFit(int nPoints, int offset, const std::vector& x f2 = roFunc(nPoints, offset, x, y, b2, aa); } sigb = .01f * sigb; - while (fabs(b2 - b1) > sigb) { + while (std::abs(b2 - b1) > sigb) { bb = b1 + .5f * (b2 - b1); if (bb == b1 || bb == b2) { break; @@ -1293,9 +1293,9 @@ float TrackResiduals::roFunc(int nPoints, int offset, const std::vector& for (int j = nPoints; j-- > 0;) { float d = y[j + offset] - (b * x[j + offset] + aa); if (y[j + offset] != 0.f) { - d /= fabs(y[j + offset]); + d /= std::abs(y[j + offset]); } - if (fabs(d) > sFloatEps) { + if (std::abs(d) > sFloatEps) { sum += (d >= 0.f ? x[j + offset] : -x[j + offset]); } } @@ -1391,7 +1391,7 @@ float TrackResiduals::getMAD2Sigma(std::vector data) const // fill vector with absolute deviations to median for (auto& entry : data) { - entry = fabs(entry - medianOfData); + entry = std::abs(entry - medianOfData); } // calculate median of abs deviations @@ -1409,7 +1409,55 @@ float TrackResiduals::getMAD2Sigma(std::vector data) const return k * medianOfAbsDeviations; } -void TrackResiduals::fitCircle(int nCl, std::array& x, std::array& y, float& xc, float& yc, float& r, std::array& residHelixY) +void TrackResiduals::fitCircle(TrackInterpolation::TrackValidationData& params) +{ + // this fast algebraic circle fit is described here: + // https://dtcenter.org/met/users/docs/write_ups/circle_fit.pdf + double xMean = 0., yMean = 0.; + int ncl = params.points.size(); + for (const auto& pnt : params.points) { + xMean += pnt.xLab; + yMean += pnt.yLab; + } + xMean /= ncl; + yMean /= ncl; + // define sums needed for circular fit + double su2 = 0., sv2 = 0., suv = 0., su3 = 0., sv3 = 0., su2v = 0., suv2 = 0.; + for (const auto& pnt : params.points) { + double ui = pnt.xLab - xMean; + double vi = pnt.yLab - yMean; + double ui2 = ui * ui; + double vi2 = vi * vi; + suv += ui * vi; + su2 += ui2; + sv2 += vi2; + su3 += ui2 * ui; + sv3 += vi2 * vi; + su2v += ui2 * vi; + suv2 += ui * vi2; + } + double rhsU = .5f * (su3 + suv2); + double rhsV = .5f * (sv3 + su2v); + double det = su2 * sv2 - suv * suv; + double uc = (rhsU * sv2 - rhsV * suv) / det; + double vc = (su2 * rhsV - suv * rhsU) / det; + double r2 = uc * uc + vc * vc + (su2 + sv2) / ncl; + params.xcLab = uc + xMean; + params.ycLab = vc + yMean; + params.r = sqrt(r2); + // write residuals to residHelixY + for (auto& pnt : params.points) { + double dx = pnt.xLab - params.xcLab; + double dxr = r2 - dx * dx; + double ys = dxr > 0 ? sqrt(dxr) : 0.f; // distance of point in y from the circle center (using fit results for r and xc) + double dy = pnt.yLab - params.ycLab; // distance of point in y from the circle center (using fit result for yc) + double dysp = dy - ys; + double dysm = dy + ys; + pnt.residHelixY = std::abs(dysp) < std::abs(dysm) ? dysp : dysm; + } +} + +void fitCircle(int nCl, std::array& x, std::array& y, float& xc, float& yc, float& r, std::array& residHelixY) { // this fast algebraic circle fit is described here: // https://dtcenter.org/met/users/docs/write_ups/circle_fit.pdf @@ -1454,11 +1502,38 @@ void TrackResiduals::fitCircle(int nCl, std::array& x, s float dy = y[i] - yc; // distance of point in y from the circle center (using fit result for yc) float dysp = dy - ys; float dysm = dy + ys; - residHelixY[i] = fabsf(dysp) < fabsf(dysm) ? dysp : dysm; + residHelixY[i] = std::abs(dysp) < std::abs(dysm) ? dysp : dysm; } // printf("r = %.4f m, xc = %.4f, yc = %.4f\n", r/100.f, xc, yc); } +bool TrackResiduals::fitPoly1(TrackInterpolation::TrackValidationData& params) +{ + // fit a straight line y = ax + b to a given set of points (x,y) + // no measurement errors assumed, no fit errors calculated + // res[0] = a (slope) + // res[1] = b (offset) + int ncl = params.points.size(); + if (ncl < 2) { + // not enough points + return false; + } + double sumX = 0., sumY = 0., sumXY = 0., sumX2 = 0., nInv = 1. / ncl; + for (const auto& pnt : params.points) { + sumX += pnt.sPath; + sumY += pnt.zTrk; + sumXY += pnt.sPath * pnt.zTrk; + sumX2 += pnt.sPath * pnt.sPath; + } + auto det = sumX2 - nInv * sumX * sumX; + if (std::abs(det) < 1e-12) { + return false; + } + params.tgl = (sumXY - nInv * sumX * sumY) / det; + params.zOffs = nInv * sumY - nInv * params.tgl * sumX; + return true; +} + bool TrackResiduals::fitPoly1(int nCl, std::array& x, std::array& y, std::array& res) { // fit a straight line y = ax + b to a given set of points (x,y) @@ -1477,7 +1552,7 @@ bool TrackResiduals::fitPoly1(int nCl, std::array& x, st sumX2 += x[i] * x[i]; } float det = sumX2 - nInv * sumX * sumX; - if (fabsf(det) < 1e-12f) { + if (std::abs(det) < 1e-12f) { return false; } res[0] = (sumXY - nInv * sumX * sumY) / det; From bc7332247248ba3bd8e455de37707b1f3e291d4b Mon Sep 17 00:00:00 2001 From: shahoian Date: Mon, 15 Jun 2026 13:43:36 +0200 Subject: [PATCH 2/2] Suppres compareToHelix and min.ref. |q/pT| cuts for extrapolated tracks By construction they were designed to eliminate bad matches of ITS/TPC to TRD and TOF --- .../include/SpacePoints/TrackInterpolation.h | 2 +- .../calibration/SpacePoints/src/TrackInterpolation.cxx | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h index 143e9679ee980..690db1b9262fb 100644 --- a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h +++ b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h @@ -353,7 +353,7 @@ class TrackInterpolation /// \param trk The track parameters, e.g. q/pT, eta, ... /// \param params Structure with per pad information recalculated on the fly /// \return 0 if the track could be validated, otherwise returns rejection code - int8_t validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes); + int8_t validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes, bool interpol); /// Filter out individual outliers from all cluster residuals of given track /// \return true for tracks which pass the cuts on e.g. max. masked clusters and false for rejected tracks diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx index abb06292fe25d..e0c3448f94629 100644 --- a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx +++ b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx @@ -701,7 +701,7 @@ void TrackInterpolation::interpolateTrack(int iSeed) mTrackValidation.clear(); // for refitted track parameters and flagging rejected clusters bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals); + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, true); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0; int iRow = 0; @@ -1043,7 +1043,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } bool stored = false; - trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals); + trackData.filterFlag = mParams->skipOutlierFiltering ? -1 : validateTrack(trackData, mTrackValidation, clusterResiduals, false); if (trackData.filterFlag <= 0 || mParams->writeUnfiltered) { int nClValidated = 0, iRow = 0; unsigned int iCl = 0; @@ -1240,7 +1240,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed) } } -int8_t TrackInterpolation::validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes) +int8_t TrackInterpolation::validateTrack(const TrackData& trk, TrackValidationData& params, const std::vector& clsRes, bool interpol) { int8_t status = 0; while (true) { @@ -1254,14 +1254,14 @@ int8_t TrackInterpolation::validateTrack(const TrackData& trk, TrackValidationDa } bool resHelix = compareToHelix(trk, params, clsRes); - if (!resHelix) { + if (!resHelix && interpol) { LOG(debug) << "Skipping track too far from helix approximation"; status |= 0x1 << 1; if (!mParams->keepRejectedResiduals) { break; // we don't keep de-validated tracks, no need to check further } } - if (std::abs(mBz) > 0.01 && std::abs(params.qpt) > mParams->maxQ2Pt) { + if (interpol && (std::abs(mBz) > 0.01 && std::abs(params.qpt) > mParams->maxQ2Pt)) { LOG(debug) << "Skipping track with too high q/pT: " << params.qpt; status |= 0x1 << 2; if (!mParams->keepRejectedResiduals) {