Skip to content

Commit 984060e

Browse files
glegrasglegras
committed
simplify error parametrization
1 parent 3c5ceb0 commit 984060e

File tree

5 files changed

+26
-49
lines changed

5 files changed

+26
-49
lines changed

‎Detectors/TRD/base/src/TrackletTransformer.cxx‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,6 @@ float TrackletTransformer::calculateDy(int detector, int slope, const PadPlane*
4040
{
4141
double padWidth = padPlane->getWidthIPad();
4242

43-
// float vDrift = mCalVdriftExB->getVdrift(detector, true);
44-
// float exb = mCalVdriftExB->getExB(detector, true);
4543
float vDrift = mVdrift[detector];
4644
float exb = mExB[detector];
4745

‎GPU/GPUTracking/DataTypes/GPUTRDRecoParam.cxx‎

Lines changed: 5 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -28,50 +28,33 @@ void GPUTRDRecoParam::init(float bz, const GPUSettingsRec* rec)
2828
if (CAMath::Abs(CAMath::Abs(bz) - 2) < 0.1) {
2929
if (bz > 0) {
3030
// magnetic field +0.2 T
31-
mRPhiA2 = resRPhiIdeal2;
32-
mRPhiB = -1.43e-2f;
3331
mRPhiC2 = 4.55e-2f;
34-
35-
mDyB = 0.035f;
36-
mCorrYDyB = 0.035f;
3732
} else {
3833
// magnetic field -0.2 T
39-
mRPhiA2 = resRPhiIdeal2;
40-
mRPhiB = 1.43e-2f;
4134
mRPhiC2 = 4.55e-2f;
42-
43-
mDyB = -0.065f;
44-
mCorrYDyB = -0.065f;
4535
}
4636
} else if (CAMath::Abs(CAMath::Abs(bz) - 5) < 0.1) {
4737
if (bz > 0) {
4838
// magnetic field +0.5 T
49-
mRPhiA2 = resRPhiIdeal2;
50-
mRPhiB = 0.125f;
5139
mRPhiC2 = 0.0961f;
52-
53-
mDyB = 0.11f;
54-
mCorrYDyB = 0.11f;
5540
} else {
5641
// magnetic field -0.5 T
57-
mRPhiA2 = resRPhiIdeal2;
58-
mRPhiB = -0.14f;
5942
mRPhiC2 = 0.1156f;
60-
61-
mDyB = -0.14f;
62-
mCorrYDyB = -0.14f;
6343
}
6444
} else {
6545
LOGP(warning, "No error parameterization available for Bz= {}. Keeping default value (sigma_y = const. = 1cm)", bz);
6646
}
47+
48+
mRPhiA2 = resRPhiIdeal2;
49+
mLorentzAngle = -0.02f + 0.13f * bz / 5.f;
6750

6851
mDyA2 = 6e-3f;
6952
mDyC2 = 0.3f;
7053
mCorrYDyA = 0.27f;
7154
mCorrYDyC = -0.44f;
7255

7356
LOGP(info, "Loaded parameterizations for Bz={}: PhiRes:[{},{},{}] DyRes:[{},{},{}] CorrYDy:[{},{},{}]",
74-
bz, mRPhiA2, mRPhiB, mRPhiC2, mDyA2, mDyB, mDyC2, mCorrYDyA, mCorrYDyB, mCorrYDyC);
57+
bz, mRPhiA2, mLorentzAngle, mRPhiC2, mDyA2, mLorentzAngle, mDyC2, mCorrYDyA, mLorentzAngle, mCorrYDyC);
7558
}
7659

7760
void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const
@@ -85,13 +68,4 @@ void GPUTRDRecoParam::recalcTrkltCov(const float tilt, const float snp, const fl
8568
cov[2] = c2 * (t2 * sy2 + sz2);
8669
}
8770

88-
void GPUTRDRecoParam::recalcTrkltCovDy(const float tilt, const float snp, float* cov) const
89-
{
90-
float t2 = tilt * tilt; // tan^2 (tilt)
91-
float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
92-
float sy2 = getRPhiRes(snp);
93-
float sdy2 = getDyRes(snp);
94-
cov[3] = getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2);
95-
cov[4] = -tilt * getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2);
96-
cov[5] = sdy2;
97-
}
71+

‎GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h‎

Lines changed: 5 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -43,42 +43,34 @@ class GPUTRDRecoParam
4343
{
4444
recalcTrkltCov(tilt, snp, rowSize, cov.data());
4545
}
46-
// covariance terms which use tracklet deflection
47-
GPUd() void recalcTrkltCovDy(const float tilt, const float snp, std::array<float, 6>& cov) const
48-
{
49-
recalcTrkltCovDy(tilt, snp, cov.data());
50-
}
5146
#endif
5247
GPUd() void recalcTrkltCov(const float tilt, const float snp, const float rowSize, float* cov) const;
53-
GPUd() void recalcTrkltCovDy(const float tilt, const float snp, float* cov) const;
5448

5549
/// Get tracklet r-phi resolution for given phi angle
5650
/// Resolution depends on the track angle sin(phi) = snp and is approximated by the formula
5751
/// sigma_y(snp) = sqrt(a^2 + c^2 * (snp - b)^2)
5852
/// more details are given in http://cds.cern.ch/record/2724259 in section 5.3.3
5953
/// \param phi angle of related track
6054
/// \return sigma_y^2 of tracklet
61-
GPUd() float getRPhiRes(float snp) const { return (mRPhiA2 + mRPhiC2 * (snp - mRPhiB) * (snp - mRPhiB)); }
62-
GPUd() float getDyRes(float snp) const { return mDyA2 + mDyC2 * (snp - mDyB) * (snp - mDyB); } // // a^2 + c^2 * (snp - b)^2
63-
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
64-
GPUd() float getCorrYDy(float snp) const { return mCorrYDyA + mCorrYDyC * (snp - mCorrYDyB) * (snp - mCorrYDyB); } // // a + c * (snp - b)^2
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
6559

6660
/// Get tracklet z correction coefficient for track-eta based corraction
6761
GPUd() float getZCorrCoeffNRC() const { return mZCorrCoefNRC; }
6862

6963
private:
7064
// tracklet error parameterization depends on the magnetic field
65+
float mLorentzAngle{0.f};
7166
// rphi
7267
float mRPhiA2{1.f}; ///< parameterization for tracklet position resolution
73-
float mRPhiB{0.f}; ///< parameterization for tracklet position resolution
7468
float mRPhiC2{0.f}; ///< parameterization for tracklet position resolution
7569
// angle
7670
float mDyA2{1.225e-3f}; ///< parameterization for tracklet angular resolution
77-
float mDyB{0.f}; ///< parameterization for tracklet angular resolution
7871
float mDyC2{0.f}; ///< parameterization for tracklet angular resolution
7972
// correlation coefficient between y residual and dy residual
8073
float mCorrYDyA{0.f};
81-
float mCorrYDyB{0.f};
8274
float mCorrYDyC{0.f};
8375

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

‎GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -595,17 +595,17 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
595595
float trkltCovTmp[3] = {0.f};
596596
if ((CAMath::Abs(deltaY) < roadY) && (CAMath::Abs(deltaZ) < roadZ)) { // TODO: check if this is still necessary after the cut before propagation of track
597597
// tracklet is in windwow: get predicted chi2 for update and store tracklet index if best guess
598-
mRecoParam->recalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmp);
598+
RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(tracklets[trkltIdx].GetZbin()), trkltCovTmp);
599599
float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp);
600600
if (Param().rec.trd.addDeflectionInChi2 && (trkWork->getSnp() < 1.f - 1e-6f) && (trkWork->getSnp() > -1.f + 1e-6f)) {
601601
// we add the slope in the chi2 calculation
602602
float trkltCovTmpWithDy[6] = {trkltCovTmp[0], trkltCovTmp[1], trkltCovTmp[2], 0.f, 0.f, 0.f};
603-
mRecoParam->recalcTrkltCovDy(tilt, trkWork->getSnp(), trkltCovTmpWithDy);
603+
RecalcTrkltCovDy(tilt, trkWork->getSnp(), trkltCovTmpWithDy);
604604
trkltCovTmpWithDy[0] += trkWork->getSigmaY2();
605605
trkltCovTmpWithDy[1] += trkWork->getSigmaZY();
606606
trkltCovTmpWithDy[2] += trkWork->getSigmaZ2();
607607

608-
// For now, dy uncertainty also includes track uncertainty, so no need to add additional uncertainty
608+
// For now, dy uncertainty parametrization also includes track uncertainty, so no need to add additional uncertainty
609609
if (InvertCov(trkltCovTmpWithDy)) {
610610
float deltaDy = spacePoints[trkltIdx].getDy() + dyTiltCorr - mRecoParam->convertAngleToDy(trkWork->getSnp());
611611
chi2 = deltaY * trkltCovTmpWithDy[0] * deltaY + 2 * deltaY * trkltCovTmpWithDy[1] * deltaZ + 2 * deltaY * trkltCovTmpWithDy[3] * deltaDy + deltaZ * trkltCovTmpWithDy[2] * deltaZ + 2 * deltaZ * trkltCovTmpWithDy[4] * deltaDy + deltaDy * trkltCovTmpWithDy[5] * deltaDy;
@@ -954,6 +954,18 @@ GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::RecalcTrkltCov(const float tilt, cons
954954
cov[2] = c2 * (t2 * sy2 + sz2);
955955
}
956956

957+
template <class TRDTRK, class PROP>
958+
GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::RecalcTrkltCovDy(const float tilt, const float snp, float (&cov)[6])
959+
{
960+
float t2 = tilt * tilt; // tan^2 (tilt)
961+
float c2 = 1.f / (1.f + t2); // cos^2 (tilt)
962+
float sy2 = mRecoParam->getRPhiRes(snp);
963+
float sdy2 = mRecoParam->getDyRes(snp);
964+
cov[3] = mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2);
965+
cov[4] = -tilt * mRecoParam->getCorrYDy(snp) * CAMath::Sqrt(sdy2 * c2 * sy2);
966+
cov[5] = sdy2;
967+
}
968+
957969
template <class TRDTRK, class PROP>
958970
GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::InvertCov(float (&cov)[6])
959971
{

‎GPU/GPUTracking/TRDTracking/GPUTRDTracker.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,7 @@ class GPUTRDTracker_t : public GPUProcessor
117117
GPUd() float GetAlphaOfSector(const int32_t sec) const;
118118
GPUd() float GetAngularPull(float dYtracklet, float snp) const;
119119
GPUd() void RecalcTrkltCov(const float tilt, const float snp, const float rowSize, float (&cov)[3]);
120+
GPUd() void RecalcTrkltCovDy(const float tilt, const float snp, float (&cov)[6]);
120121
GPUd() bool InvertCov(float (&cov)[6]);
121122
GPUd() void FindChambersInRoad(const TRDTRK* t, const float roadY, const float roadZ, const int32_t iLayer, int32_t* det, const float zMax, const float alpha, const float zShiftTrk) const;
122123
GPUd() bool IsGeoFindable(const TRDTRK* t, const int32_t layer, const float alpha, const float zShiftTrk) const;

0 commit comments

Comments
 (0)