Skip to content

Commit 501f33d

Browse files
committed
update error parametrization
1 parent db38fa1 commit 501f33d

7 files changed

Lines changed: 101 additions & 49 deletions

File tree

Detectors/TRD/qc/src/Tracking.cxx

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,7 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
149149
const PadPlane* pad = Geometry::instance()->getPadPlane(trkltDet);
150150
float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees
151151
float tiltCorrUp = tilt * (mTrackletsCalib[trkltId].getZ() - trk.getZ());
152+
float dyTiltCorr = tilt * trk.getTgl() * Geometry::instance()->cdrHght();
152153
float zPosCorrUp = mTrackletsCalib[trkltId].getZ() + mRecoParam.getZCorrCoeffNRC() * trk.getTgl();
153154
float padLength = pad->getRowSize(tracklet.getPadRow());
154155
if (!((trk.getSigmaZ2() < (padLength * padLength / 12.f)) && (std::fabs(mTrackletsCalib[trkltId].getZ() - trk.getZ()) < padLength))) {
@@ -159,10 +160,12 @@ void Tracking::checkTrack(const TrackTRD& trkTrd, bool isTPCTRD)
159160
float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f;
160161
float yCorrPileUp = tCorrPileUp * slopeFactor;
161162
float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
163+
164+
float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trk.getSnp())) / std::sqrt(mRecoParam.getDyRes(trk.getSnp(), 0));
162165

163166
std::array<float, 2> trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp};
164167
std::array<float, 3> trkltCovUp;
165-
mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp);
168+
mRecoParam.recalcTrkltCov(tilt, trk.getSnp(), pad->getRowSize(tracklet.getPadRow()), trkltCovUp, angularPull, 0);
166169
trkltCovUp[0] += yAddErrPileUp2;
167170
auto chi2trklt = trk.getPredictedChi2(trkltPosUp, trkltCovUp);
168171

Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -879,6 +879,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
879879
}
880880
const PadPlane* pad = Geometry::instance()->getPadPlane(trkltDet);
881881
float tilt = tan(TMath::DegToRad() * pad->getTiltingAngle()); // tilt is signed! and returned in degrees
882+
float dyTiltCorr = tilt * trkParam->getTgl() * Geometry::instance()->cdrHght();
882883
float tiltCorrUp = tilt * (mTrackletsCalib[trkltId].getZ() - trkParam->getZ());
883884
float zPosCorrUp = mTrackletsCalib[trkltId].getZ() + mRecoParam.getZCorrCoeffNRC() * trkParam->getTgl();
884885
float padLength = pad->getRowSize(mTrackletsRaw[trkltId].getPadRow());
@@ -890,10 +891,13 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,
890891
float slopeFactor = mTrackletsRaw[trkltId].getSlopeFloat() * pad->getWidthIPad() / 4.f;
891892
float yCorrPileUp = tCorrPileUp * slopeFactor;
892893
float yAddErrPileUp2 = tErrPileUp2 * slopeFactor * slopeFactor;
894+
895+
int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet);
896+
float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber));
893897

894898
std::array<float, 2> trkltPosUp{mTrackletsCalib[trkltId].getY() - tiltCorrUp + yCorrPileUp, zPosCorrUp};
895899
std::array<float, 3> trkltCovUp;
896-
mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp);
900+
mRecoParam.recalcTrkltCov(tilt, trkParam->getSnp(), pad->getRowSize(mTrackletsRaw[trkltId].getPadRow()), trkltCovUp, (mRec->GetParam().rec.trd.useAngularPull != 0 ? angularPull : 0.), nTrackletsChamber);
897901
trkltCovUp[0] += yAddErrPileUp2;
898902

899903
chi2 += trkParam->getPredictedChi2(trkltPosUp, trkltCovUp);

GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -24,16 +24,19 @@ using namespace o2::gpu;
2424
// error parameterizations taken from http://cds.cern.ch/record/2724259 Appendix A
2525
void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec)
2626
{
27-
float resRPhiIdeal2 = 1.6e-3f;
27+
float resRPhiIdeal = 0.04f;
28+
float resVsTanPhiMisalign = 0.f;
2829
if (rec) {
29-
resRPhiIdeal2 = rec->trd.trkltResRPhiIdeal * rec->trd.trkltResRPhiIdeal;
30+
resRPhiIdeal = rec->trd.trkltResRPhiIdeal;
31+
resVsTanPhiMisalign = rec->trd.trkltResVsTanPhiMisalign;
3032
mPileUpRangeBefore = -rec->trd.pileupBwdNBC;
3133
mPileUpRangeAfter = rec->trd.pileupFwdNBC;
3234
}
3335
#ifndef GPUCA_STANDALONE
3436
else {
3537
const auto& rtrd = GPU_GET_CONFIG(GPUSettingsRecTRD);
36-
resRPhiIdeal2 = rtrd.trkltResRPhiIdeal * rtrd.trkltResRPhiIdeal;
38+
resRPhiIdeal = rtrd.trkltResRPhiIdeal;
39+
resVsTanPhiMisalign = rtrd.trkltResVsTanPhiMisalign;
3740
mPileUpRangeBefore = -rtrd.pileupBwdNBC;
3841
mPileUpRangeAfter = rtrd.pileupFwdNBC;
3942
}
@@ -59,23 +62,22 @@ void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec)
5962
LOGP(warning, "No error parameterization available for Bz= {}. Keeping default value (sigma_y = const. = 1cm)", bz);
6063
}
6164

62-
mRPhiA2 = resRPhiIdeal2;
65+
mRPhiA = resRPhiIdeal;
66+
mRPhiATgp = resVsTanPhiMisalign;
6367
mLorentzAngle = -0.02f + 0.13f * bz / 5.f;
6468

6569
mDyA2 = 6e-3f;
6670
mDyC2 = 0.3f;
67-
mCorrYDyA = 0.27f;
68-
mCorrYDyC = -0.44f;
6971

70-
LOGP(info, "Loaded parameterizations for Bz={}: PhiRes:[{},{},{}] DyRes:[{},{},{}] CorrYDy:[{},{},{}]",
71-
bz, mRPhiA2, mLorentzAngle, mRPhiC2, mDyA2, mLorentzAngle, mDyC2, mCorrYDyA, mLorentzAngle, mCorrYDyC);
72+
LOGP(info, "Loaded parameterizations for Bz={}: PhiRes:[{},{},{},{}] DyRes:[{},{},{}]",
73+
bz, mRPhiA, mRPhiATgp, mLorentzAngle, mRPhiC2, mDyA2, mLorentzAngle, mDyC2);
7274
}
7375

74-
void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const
76+
void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov, const float pull, const int occupancy) const
7577
{
7678
float t2 = tilt * tilt; // tan^2 (tilt)
7779
float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
78-
float sy2 = getRPhiRes(snp);
80+
float sy2 = getRPhiRes(snp, CAMath::Abs(pull), occupancy);
7981
float sz2 = rowSize * rowSize / 12.f;
8082
cov[0] = c2 * (sy2 + t2 * sz2);
8183
cov[1] = c2 * tilt * (sz2 - sy2);

GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h

Lines changed: 40 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -39,23 +39,17 @@ class GPUTRDRecoParam
3939

4040
#if !defined(GPUCA_GPUCODE_DEVICE)
4141
/// Recalculate tracklet covariance based on phi angle of related track
42-
GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, std::array<float, 3>& cov) const
42+
GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, std::array<float, 3>& cov, const float pull = 0., const int occupancy = 0) const
4343
{
44-
recalcTrkltCov(tilt, snp, rowSize, cov.data());
44+
recalcTrkltCov(tilt, snp, rowSize, cov.data(), pull, occupancy);
4545
}
4646
#endif
47-
GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const;
48-
49-
/// Get tracklet r-phi resolution for given phi angle
50-
/// Resolution depends on the track angle sin(phi) = snp and is approximated by the formula
51-
/// sigma_y(snp) = sqrt(a^2 + c^2 * (snp - b)^2)
52-
/// more details are given in http://cds.cern.ch/record/2724259 in section 5.3.3
53-
/// \param phi angle of related track
54-
/// \return sigma_y^2 of tracklet
55-
GPUd() float getRPhiRes(float snp) const { return (mRPhiA2 + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle)); }
56-
GPUd() float getDyRes(float snp) const { return mDyA2 + mDyC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle); } // a^2 + c^2 * (snp - b)^2
57-
GPUd() float convertAngleToDy(float snp) const { return 3.f * snp / CAMath::Sqrt(1 - snp * snp); } // when calibrated, sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) works well
58-
GPUd() float getCorrYDy(float snp) const { return mCorrYDyA + mCorrYDyC * (snp - mLorentzAngle) * (snp - mLorentzAngle); } // a + c * (snp - b)^2
47+
GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov, const float pull = 0., const int occupancy = 0) const;
48+
49+
GPUd() float getRPhiRes(float snp, float pull = 0.f, int occupancy = 0) const;
50+
GPUd() float getDyRes(float snp, int occupancy = 0) const { return mDyA2 + mDyC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + mOccDyA * occupancy; } // a^2 + c^2 * (snp - b)^2
51+
GPUd() float convertAngleToDy(float snp) const { return 3.f * snp / CAMath::Sqrt(1 - snp * snp); } // when calibrated, sin(phi) = (dy / xDrift) / sqrt(1+(dy/xDrift)^2) works well
52+
GPUd() float getCorrYDy() const { return mCorrYDy; }
5953
GPUd() float getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0 = true, bool Q1 = true) const;
6054
GPUd() float getPileUpProbTrack(int nBC, std::array<int, 6> Q0, std::array<int, 6> Q1) const;
6155

@@ -70,14 +64,21 @@ class GPUTRDRecoParam
7064
// tracklet error parameterization depends on the magnetic field
7165
float mLorentzAngle{0.f};
7266
// rphi
73-
float mRPhiA2{1.f}; ///< parameterization for tracklet position resolution
74-
float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution
67+
float mRPhiA{1.f}; ///< parameterization for tracklet position resolution
68+
float mRPhiATgp{1.f}; ///< parameterization for tracklet position resolution
69+
float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution
7570
// angle
7671
float mDyA2{1.225e-3f}; ///< parameterization for tracklet angular resolution
7772
float mDyC2{0.f}; ///< parameterization for tracklet angular resolution
78-
// correlation coefficient between y residual and dy residual
79-
float mCorrYDyA{0.f};
80-
float mCorrYDyC{0.f};
73+
// variation in y when dy variates by one sigma (= cov / sigma_dy = corr * sigma_y) (valid within 2sigma of dy)
74+
float mCorrYDy{0.13f};
75+
// error parametrization vs angular pull (pol2)
76+
float mPullA{6.8e-3f};
77+
float mPullB{0.049f};
78+
// error parametrization of y position vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy))
79+
float mOccA{3.3e-4f};
80+
// error parametrization for dy vs occupancy defined as ntracklets within chamber (prop to sqrt(occupancy))
81+
float mOccDyA{2.5e-4f};
8182

8283
float mZCorrCoefNRC{1.4f}; ///< tracklet z-position depends linearly on track dip angle
8384

@@ -106,6 +107,26 @@ class GPUTRDRecoParam
106107
ClassDefNV(GPUTRDRecoParam, 4);
107108
};
108109

110+
/// Get tracklet r-phi resolution for given phi angle
111+
/// Resolution depends on the track angle sin(phi) = snp and is approximated by the formula
112+
/// sigma_y(snp) = sqrt(a^2 + c^2 * (snp - b)^2)
113+
/// more details are given in http://cds.cern.ch/record/2724259 in section 5.3.3
114+
/// \param phi angle of related track
115+
/// \return sigma_y^2 of tracklet
116+
/// also depend on absolute pull and on chamber occupancy
117+
GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy) const {
118+
// flat uncertainty + radial-alignment uncertainty depending on tan(phi)
119+
float tgp = (CAMath::Abs(snp) < 0.99999f) ? CAMath::Abs(snp) / CAMath::Sqrt(1 - snp * snp) : 1e6;
120+
float resIdeal = mRPhiA + mRPhiATgp * tgp;
121+
if (pull > 10) {
122+
// parametrization does not really work well for such large pull values
123+
pull = 10.f;
124+
}
125+
float resPull = mPullA * pull * pull + mPullB * pull; // parametrization as pol2 summed in quadrature
126+
float resOccupancy = mOccA * occupancy; // parametrization as sqrt() summed in quadrature
127+
return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy);
128+
}
129+
109130
GPUdi() float GPUTRDRecoParam::getPileUpProbTracklet(int nBC, bool withChargeInfo, bool Q0, bool Q1) const
110131
{
111132
// get the probability that the tracklet with charges Q0 and Q1 belongs to a given BC, with a (signed) distance nBC from the TRD-triggered BC

GPU/GPUTracking/Definitions/GPUSettingsList.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -182,6 +182,7 @@ AddOptionRTC(addTimeRoadITSTPC, float, 2.5f, "", 0, "Increase time search road b
182182
AddOptionRTC(extraRoadY, float, 5.f, "", 0, "Addition to search road around track prolongation along Y in cm")
183183
AddOptionRTC(extraRoadZ, float, 10.f, "", 0, "Addition to search road around track prolongation along Z in cm")
184184
AddOptionRTC(trkltResRPhiIdeal, float, 1.f, "", 0, "Optimal tracklet rphi resolution in cm (in case phi of track = lorentz angle)")
185+
AddOptionRTC(trkltResVsTanPhiMisalign, float, 0.f, "", 0, "tan(phi) dependence of tracklet error due to radial misalignment (centered at 0 angle)")
185186
AddOptionRTC(maxChi2Red, float, 99.f, "", 0, "maximum chi2 per attached tracklet for TRD tracks TODO: currently effectively disabled, requires tuning")
186187
AddOptionRTC(applyDeflectionCut, uint8_t, 0, "", 0, "Set to 1 to enable tracklet selection based on deflection")
187188
AddOptionRTC(addDeflectionInChi2, uint8_t, 0, "", 0, "Set to 1 to add the deflection in the chi2 calculation for matching")
@@ -190,6 +191,7 @@ AddOptionRTC(nTrackletsMin, uint8_t, 3, "", 0, "Tracks with less attached trackl
190191
AddOptionRTC(matCorrType, uint8_t, 2, "", 0, "Material correction to use: 0 - none, 1 - TGeo, 2 - matLUT")
191192
AddOptionRTC(pileupFwdNBC, uint8_t, 70, "", 0, "Post-trigger Pile-up integration time in BCs")
192193
AddOptionRTC(pileupBwdNBC, uint8_t, 130, "", 0, "Pre-trigger Pile-up integration time in BCs")
194+
AddOptionRTC(useAngularPull, uint8_t, 1, "", 0, "0 = don't use angular pull; 1 = additional error based on angular pull for refit only; 2 = add error also for chi2")
193195
AddHelp("help", 'h')
194196
EndConfig()
195197

0 commit comments

Comments
 (0)