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, 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 - 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..e0c3448f94629 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, true); + 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, false); + 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, bool interpol) { - 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 && 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 (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) { + 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;