Skip to content

Commit a866387

Browse files
wiechulasawenzel
wiechula
authored and
sawenzel
committed
Move check of outliers in residuals to proper place
1 parent b0f13b6 commit a866387

File tree

3 files changed

+76
-61
lines changed

3 files changed

+76
-61
lines changed

‎Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibParam.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ static constexpr int NZ2XBins = 5; ///< number of bins in z/x
5656
static constexpr float MaxResid = 20.f; ///< max residual in y and z
5757
static constexpr float MaxY = 50.f; ///< max value for y position (sector coordinates)
5858
static constexpr float MaxZ = 300.f; ///< max value for z position
59-
static constexpr float MaxTgSlp = 1.f; ///< max value for phi (from snp)
59+
static constexpr float MaxTgSlp = 1.f; ///< max value for phi (from snp, converted to tangens)
6060

6161
// miscellaneous
6262
static constexpr float sEps = 1e-6f; ///< small number for float comparisons

‎Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h‎

Lines changed: 30 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -60,11 +60,6 @@ struct TPCClusterResiduals {
6060
float snp{}; ///< sin of the phi angle between padrow and track
6161
unsigned char sec{}; ///< sector number 0..35
6262
unsigned char dRow{}; ///< distance to previous row in units of pad rows
63-
void setDY(float val) { dy = fabs(val) < param::MaxResid ? val : std::copysign(param::MaxResid, val); }
64-
void setDZ(float val) { dz = fabs(val) < param::MaxResid ? val : std::copysign(param::MaxResid, val); }
65-
void setY(float val) { y = fabs(val) < param::MaxY ? val : std::copysign(param::MaxY, val); }
66-
void setZ(float val) { z = fabs(val) < param::MaxZ ? val : std::copysign(param::MaxZ, val); }
67-
void setSnp(float val) { snp = fabs(val) < param::MaxTgSlp ? val : std::copysign(param::MaxTgSlp, val); }
6863
ClassDefNV(TPCClusterResiduals, 4);
6964
};
7065

@@ -120,16 +115,16 @@ struct TrackDataExtended {
120115

121116
/// Structure filled for each track with track quality information and a vector with TPCClusterResiduals
122117
struct TrackData {
123-
o2::dataformats::GlobalTrackID gid{}; ///< global track ID for seeding track
124-
o2::track::TrackPar par{}; ///< ITS track at inner TPC radius
125-
float dEdxTPC{}; ///< TPC dEdx information
126-
float chi2TPC{}; ///< chi2 of TPC track
127-
float chi2ITS{}; ///< chi2 of ITS track
128-
float chi2TRD{}; ///< chi2 of TRD track
129-
unsigned short nClsTPC{}; ///< number of attached TPC clusters
130-
unsigned short nClsITS{}; ///< number of attached ITS clusters
131-
unsigned short nTrkltsTRD{}; ///< number of attached TRD tracklets
132-
unsigned short clAvailTOF{}; ///< whether or not track seed has a matched TOF cluster
118+
o2::dataformats::GlobalTrackID gid{}; ///< global track ID for seeding track
119+
o2::track::TrackPar par{}; ///< ITS track at inner TPC radius
120+
float dEdxTPC{}; ///< TPC dEdx information
121+
float chi2TPC{}; ///< chi2 of TPC track
122+
float chi2ITS{}; ///< chi2 of ITS track
123+
float chi2TRD{}; ///< chi2 of TRD track
124+
unsigned short nClsTPC{}; ///< number of attached TPC clusters
125+
unsigned short nClsITS{}; ///< number of attached ITS clusters
126+
unsigned short nTrkltsTRD{}; ///< number of attached TRD tracklets
127+
unsigned short clAvailTOF{}; ///< whether or not track seed has a matched TOF cluster
133128
o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track
134129
ClassDefNV(TrackData, 6);
135130
};
@@ -282,30 +277,30 @@ class TrackInterpolation
282277
static constexpr float sFloatEps{1.e-7f}; ///< float epsilon for robust linear fitting
283278
// parameters + settings
284279
const SpacePointsCalibConfParam* mParams = nullptr;
285-
float mTPCTimeBinMUS{.2f}; ///< TPC time bin duration in us
286-
float mTPCVDriftRef = -1.; ///< TPC nominal drift speed in cm/microseconds
287-
float mTPCDriftTimeOffsetRef = 0.; ///< TPC nominal (e.g. at the start of run) drift time bias in cm/mus
288-
float mSqrtS{13600.f}; ///< centre of mass energy set from LHC IF
289-
MatCorrType mMatCorr{MatCorrType::USEMatCorrNONE}; ///< if material correction should be done
290-
int mMaxTracksPerTF{-1}; ///< max number of tracks to be processed per TF (-1 means there is no limit)
291-
int mAddTracksForMapPerTF{0}; ///< in case residuals from different track types are used for vDrift calibration and map creation this defines the statistics for the latter
292-
bool mDumpTrackPoints{false}; ///< dump also track points in ITS, TRD and TOF
293-
bool mProcessSeeds{false}; ///< in case for global tracks also their shorter parts are processed separately
294-
bool mProcessITSTPConly{false}; ///< flag, whether or not to extrapolate ITS-only through TPC
280+
float mTPCTimeBinMUS{.2f}; ///< TPC time bin duration in us
281+
float mTPCVDriftRef = -1.; ///< TPC nominal drift speed in cm/microseconds
282+
float mTPCDriftTimeOffsetRef = 0.; ///< TPC nominal (e.g. at the start of run) drift time bias in cm/mus
283+
float mSqrtS{13600.f}; ///< centre of mass energy set from LHC IF
284+
MatCorrType mMatCorr{MatCorrType::USEMatCorrNONE}; ///< if material correction should be done
285+
int mMaxTracksPerTF{-1}; ///< max number of tracks to be processed per TF (-1 means there is no limit)
286+
int mAddTracksForMapPerTF{0}; ///< in case residuals from different track types are used for vDrift calibration and map creation this defines the statistics for the latter
287+
bool mDumpTrackPoints{false}; ///< dump also track points in ITS, TRD and TOF
288+
bool mProcessSeeds{false}; ///< in case for global tracks also their shorter parts are processed separately
289+
bool mProcessITSTPConly{false}; ///< flag, whether or not to extrapolate ITS-only through TPC
295290
o2::dataformats::GlobalTrackID::mask_t mSourcesConfigured; ///< the track sources taken into account for extra-/interpolation
296291
o2::dataformats::GlobalTrackID::mask_t mSourcesConfiguredMap; ///< possible subset of mSourcesConfigured
297292
bool mSingleSourcesConfigured{true}; ///< whether mSourcesConfigured == mSourcesConfiguredMap
298293

299294
// input
300-
const o2::globaltracking::RecoContainer* mRecoCont = nullptr; ///< input reco container
301-
std::vector<o2::dataformats::GlobalTrackID> mGIDs{}; ///< GIDs of input tracks
302-
std::vector<o2::globaltracking::RecoContainer::GlobalIDSet> mGIDtables{}; ///< GIDs of contributors from single detectors for each seed
303-
std::vector<float> mTrackTimes{}; ///< time estimates for all input tracks in micro seconds
304-
std::vector<o2::track::TrackParCov> mSeeds{}; ///< seeding track parameters (ITS tracks)
305-
std::map<int, int> mTrackTypes; ///< mapping of track source to array index in mTrackIndices
306-
std::array<std::vector<uint32_t>, 4> mTrackIndices; ///< keep GIDs of input tracks separately for each track type
307-
gsl::span<const TPCClRefElem> mTPCTracksClusIdx; ///< input TPC cluster indices from span
308-
const ClusterNativeAccess* mTPCClusterIdxStruct = nullptr; ///< struct holding the TPC cluster indices
295+
const o2::globaltracking::RecoContainer* mRecoCont = nullptr; ///< input reco container
296+
std::vector<o2::dataformats::GlobalTrackID> mGIDs{}; ///< GIDs of input tracks
297+
std::vector<o2::globaltracking::RecoContainer::GlobalIDSet> mGIDtables{}; ///< GIDs of contributors from single detectors for each seed
298+
std::vector<float> mTrackTimes{}; ///< time estimates for all input tracks in micro seconds
299+
std::vector<o2::track::TrackParCov> mSeeds{}; ///< seeding track parameters (ITS tracks)
300+
std::map<int, int> mTrackTypes; ///< mapping of track source to array index in mTrackIndices
301+
std::array<std::vector<uint32_t>, 4> mTrackIndices; ///< keep GIDs of input tracks separately for each track type
302+
gsl::span<const TPCClRefElem> mTPCTracksClusIdx; ///< input TPC cluster indices from span
303+
const ClusterNativeAccess* mTPCClusterIdxStruct = nullptr; ///< struct holding the TPC cluster indices
309304
// ITS specific input only needed for debugging
310305
gsl::span<const int> mITSTrackClusIdx; ///< input ITS track cluster indices span
311306
std::vector<o2::BaseCluster<float>> mITSClustersArray; ///< ITS clusters created in run() method from compact clusters
@@ -329,6 +324,7 @@ class TrackInterpolation
329324
std::unique_ptr<TPCFastTransform> mFastTransform{}; ///< TPC cluster transformation
330325
float mBz; ///< required for helix approximation
331326
bool mInitDone{false}; ///< initialization done flag
327+
size_t mRejectedResiduals{}; ///< number of rejected residuals
332328

333329
ClassDefNV(TrackInterpolation, 1);
334330
};

‎Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx‎

Lines changed: 45 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -334,7 +334,8 @@ void TrackInterpolation::process()
334334
extrapolateTrack(iSeed);
335335
}
336336
}
337-
LOG(info) << "Could process " << mTrackData.size() << " tracks successfully";
337+
LOG(info) << "Could process " << mTrackData.size() << " tracks successfully. " << mRejectedResiduals << " residuals were rejected. " << mClRes.size() << " residuals were accepted.";
338+
mRejectedResiduals = 0;
338339
}
339340

340341
void TrackInterpolation::interpolateTrack(int iSeed)
@@ -404,7 +405,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
404405
mCache[iRow].szy[ExtOut] = trkWork.getSigmaZY();
405406
mCache[iRow].sz2[ExtOut] = trkWork.getSigmaZ2();
406407
mCache[iRow].snp[ExtOut] = trkWork.getSnp();
407-
//printf("Track alpha at row %i: %.2f, Y(%.2f), Z(%.2f)\n", iRow, trkWork.getAlpha(), trkWork.getY(), trkWork.getZ());
408+
// printf("Track alpha at row %i: %.2f, Y(%.2f), Z(%.2f)\n", iRow, trkWork.getAlpha(), trkWork.getY(), trkWork.getZ());
408409
}
409410

410411
// start from outermost cluster with outer refit and back propagation
@@ -431,7 +432,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
431432
// TODO: check if reset of covariance matrix is needed here (or, in case TOF point is not available at outermost TRD layer)
432433
if (!trkWork.update(clTOFYZ, clTOFCov)) {
433434
LOG(debug) << "Failed to update extrapolated ITS track with TOF cluster";
434-
//LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]);
435+
// LOGF(info, "trkWork.y=%f, cl.y=%f, trkWork.z=%f, cl.z=%f", trkWork.getY(), clTOFYZ[0], trkWork.getZ(), clTOFYZ[1]);
435436
return;
436437
}
437438
}
@@ -509,7 +510,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
509510
}
510511
if (!propagator->PropagateToXBxByBz(trkWork, param::RowX[iRow], mParams->maxSnp, mParams->maxStep, mMatCorr)) {
511512
LOG(debug) << "Failed on back propagation";
512-
//printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha);
513+
// printf("trkX(%.2f), clX(%.2f), clY(%.2f), clZ(%.2f), alphaTOF(%.2f)\n", trkWork.getX(), param::RowX[iRow], clTOFYZ[0], clTOFYZ[1], clTOFAlpha);
513514
return;
514515
}
515516
mCache[iRow].y[ExtIn] = trkWork.getY();
@@ -535,15 +536,14 @@ void TrackInterpolation::interpolateTrack(int iSeed)
535536
// simple average w/o weighting for angle
536537
mCache[iRow].snp[Int] = (mCache[iRow].snp[ExtOut] + mCache[iRow].snp[ExtIn]) / 2.f;
537538

538-
TPCClusterResiduals res;
539-
res.setDY(mCache[iRow].clY - mCache[iRow].y[Int]);
540-
res.setDZ(mCache[iRow].clZ - mCache[iRow].z[Int]);
541-
res.setY(mCache[iRow].y[Int]);
542-
res.setZ(mCache[iRow].z[Int]);
543-
res.setSnp(mCache[iRow].snp[Int]);
544-
res.sec = mCache[iRow].clSec;
545-
res.dRow = deltaRow;
546-
clusterResiduals.push_back(std::move(res));
539+
const auto dY = mCache[iRow].clY - mCache[iRow].y[Int];
540+
const auto dZ = mCache[iRow].clZ - mCache[iRow].z[Int];
541+
const auto y = mCache[iRow].y[Int];
542+
const auto z = mCache[iRow].z[Int];
543+
const auto snp = mCache[iRow].snp[Int];
544+
const auto sec = mCache[iRow].clSec;
545+
clusterResiduals.emplace_back(dY, dZ, y, z, snp, sec, deltaRow);
546+
547547
deltaRow = 1;
548548
}
549549
trackData.chi2TRD = gidTable[GTrackID::TRD].isIndexSet() ? mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]).getChi2() : 0;
@@ -567,8 +567,17 @@ void TrackInterpolation::interpolateTrack(int iSeed)
567567
continue;
568568
}
569569
++nClValidated;
570-
float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
571-
mClRes.emplace_back(clusterResiduals[iCl].dy, clusterResiduals[iCl].dz, tgPhi, clusterResiduals[iCl].y, clusterResiduals[iCl].z, iRow, clusterResiduals[iCl].sec);
570+
const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
571+
const auto dy = clusterResiduals[iCl].dy;
572+
const auto dz = clusterResiduals[iCl].dz;
573+
const auto y = clusterResiduals[iCl].y;
574+
const auto z = clusterResiduals[iCl].z;
575+
const auto sec = clusterResiduals[iCl].sec;
576+
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)) {
577+
mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, sec);
578+
} else {
579+
++mRejectedResiduals;
580+
}
572581
}
573582
trackData.clIdx.setEntries(nClValidated);
574583
mTrackData.push_back(std::move(trackData));
@@ -645,16 +654,17 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
645654
if (!propagator->PropagateToXBxByBz(trkWork, x, mParams->maxSnp, mParams->maxStep, mMatCorr)) {
646655
return;
647656
}
648-
TPCClusterResiduals res;
649-
res.setDY(y - trkWork.getY());
650-
res.setDZ(z - trkWork.getZ());
651-
res.setY(trkWork.getY());
652-
res.setZ(trkWork.getZ());
653-
res.setSnp(trkWork.getSnp());
654-
res.sec = sector;
655-
res.dRow = row - rowPrev;
657+
658+
const auto dY = y - trkWork.getY();
659+
const auto dZ = z - trkWork.getZ();
660+
const auto ty = trkWork.getY();
661+
const auto tz = trkWork.getZ();
662+
const auto snp = trkWork.getSnp();
663+
const auto sec = sector;
664+
665+
clusterResiduals.emplace_back(dY, dZ, ty, tz, snp, sec, row - rowPrev);
666+
656667
rowPrev = row;
657-
clusterResiduals.push_back(std::move(res));
658668
++nMeasurements;
659669
}
660670
trackData.chi2TPC = trkTPC.getChi2();
@@ -683,8 +693,17 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
683693
continue;
684694
}
685695
++nClValidated;
686-
float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
687-
mClRes.emplace_back(clusterResiduals[iCl].dy, clusterResiduals[iCl].dz, tgPhi, clusterResiduals[iCl].y, clusterResiduals[iCl].z, iRow, clusterResiduals[iCl].sec);
696+
const float tgPhi = clusterResiduals[iCl].snp / std::sqrt((1.f - clusterResiduals[iCl].snp) * (1.f + clusterResiduals[iCl].snp));
697+
const auto dy = clusterResiduals[iCl].dy;
698+
const auto dz = clusterResiduals[iCl].dz;
699+
const auto y = clusterResiduals[iCl].y;
700+
const auto z = clusterResiduals[iCl].z;
701+
const auto sec = clusterResiduals[iCl].sec;
702+
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)) {
703+
mClRes.emplace_back(dy, dz, tgPhi, y, z, iRow, sec);
704+
} else {
705+
++mRejectedResiduals;
706+
}
688707
}
689708
trackData.clIdx.setEntries(nClValidated);
690709
mTrackData.push_back(std::move(trackData));

0 commit comments

Comments
 (0)