Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 11 additions & 13 deletions Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -153,9 +153,9 @@ void TRDGlobalTracking::updateTimeDependentParams(ProcessingContext& pc)
mBase->init(pc);
mBase->setLocalGainFactors(pc.inputs().get<o2::trd::LocalGainFactor*>("localgainfactors").get());
}

pc.inputs().get<std::array<int, constants::MAXCHAMBER>*>("chamberstatus"); // called to trigger finaliseCCDB
pc.inputs().get<o2::trd::PadStatus*>("padstatus"); // called to trigger finaliseCCDB
pc.inputs().get<o2::trd::PadStatus*>("padstatus"); // called to trigger finaliseCCDB
}

const auto& trackTune = TrackTuneParams::Instance();
Expand Down Expand Up @@ -207,9 +207,8 @@ void TRDGlobalTracking::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
if ((*chamberStatus)[iDet] == 3) {
mTracker->SetChamberStatus(iDet, false); // chamber is good
}
else {
mTracker->SetChamberStatus(iDet, true); // chamber is bad
} else {
mTracker->SetChamberStatus(iDet, true); // chamber is bad
}
}
return;
Expand All @@ -219,14 +218,13 @@ void TRDGlobalTracking::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
const o2::trd::PadStatus* padStatus = (const o2::trd::PadStatus*)obj;
for (int iDet = 0; iDet < constants::MAXCHAMBER; iDet++) {
for (int iCol = 0; iCol < constants::NCOLUMN; iCol++) {
for (int iRow = 0; iRow < ((iDet%30)/6 == 2 ? constants::NROWC0 : constants::NROWC1); iRow++) {
for (int iRow = 0; iRow < ((iDet % 30) / 6 == 2 ? constants::NROWC0 : constants::NROWC1); iRow++) {
if (padStatus->isMasked(iDet, iCol, iRow) || padStatus->isNotConnected(iDet, iCol, iRow)) {
mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, true); // pad is masked
}
else {
} else {
mTracker->SetPadStatus(iDet * constants::NCOLUMN * constants::NROWC1 + iCol * constants::NROWC1 + iRow, false); // pad is not masked
}
}
}
}
}
return;
Expand Down Expand Up @@ -533,7 +531,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
if (trdTrack.getChi2() / trdTrack.getNtracklets() > mTracker->Param().rec.trd.maxChi2Red) {
continue;
}

// Find most probable BCs and RMS for pile-up correction and error. Same BC is assumed for all tracklets
float maxProb = 0.f;
// The uncertainty is the RMS wrt the default correction of all possible corrections weighted by their probability
Expand Down Expand Up @@ -806,7 +804,7 @@ bool TRDGlobalTracking::refitTPCTRDTrack(TrackTRD& trk, float timeTRD, o2::globa
return false;
}
if (pileUpOn) { // account pileup time uncertainty in Z errors
//timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS();
// timeZErr = mTPCVdrift * trk.getPileUpTimeErrorMUS();
timeZErr = mTPCVdrift * mTPCVdrift * mTErrPileUp2;
outerParam.updateCov(timeZErr, o2::track::CovLabels::kSigZ2);
}
Expand Down Expand Up @@ -930,7 +928,7 @@ bool TRDGlobalTracking::refitTRDTrack(TrackTRD& trk, float& chi2, bool inwards,

int nTrackletsChamber = mTracker->GetNtrackletsChamber(trk.getCollisionId(), trkltDet);
float angularPull = (mTrackletsCalib[trkltId].getDy() + dyTiltCorr - mRecoParam.convertAngleToDy(trkParam->getSnp())) / std::sqrt(mRecoParam.getDyRes(trkParam->getSnp(), nTrackletsChamber));

// Correction of y position based on angular pull
if (mRec->GetParam().rec.trd.useAngularPull == 3 || mRec->GetParam().rec.trd.useAngularPull == 4) {
float corrPull = -angularPull * mRecoParam.getCorrYDy(trkParam->getSnp());
Expand Down Expand Up @@ -1046,7 +1044,7 @@ DataProcessorSpec getTRDGlobalTrackingSpec(bool useMC, GTrackID::mask_t src, boo
// request calibration data
inputs.emplace_back("localgainfactors", "TRD", "LOCALGAINFACTORS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/LocalGainFactor"));
}

// request list of bad chambers and masked pads to estimate better the number of findable tracklets
inputs.emplace_back("chamberstatus", "TRD", "CHAMBERSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/DCSDPsFedChamberStatus"));
inputs.emplace_back("padstatus", "TRD", "PADSTATUS", 0, Lifetime::Condition, ccdbParamSpec("TRD/Calib/PadStatus"));
Expand Down
17 changes: 8 additions & 9 deletions GPU/GPUTracking/DataTypes/GPUTRDRecoParam.h
Original file line number Diff line number Diff line change
Expand Up @@ -136,32 +136,31 @@ GPUdi() float GPUTRDRecoParam::getRPhiRes(float snp, float pull, int occupancy)
return (resIdeal * resIdeal + mRPhiC2 * (snp - mLorentzAngle) * (snp - mLorentzAngle) + resPull * resPull + resOccupancy);
}

GPUdi() double GPUTRDRecoParam::getDyLikelihood(float snp, float slope, int occupancy) const
{
GPUdi() double GPUTRDRecoParam::getDyLikelihood(float snp, float slope, int occupancy) const
{
// Gaussian + left exponential + right exponential
double likelihood = CAMath::Exp(- 0.5f * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) / GPUTRDRecoParam::getDyRes(snp, occupancy));
double likelihood = CAMath::Exp(-0.5f * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)) / GPUTRDRecoParam::getDyRes(snp, occupancy));

// Normalization for the exponential is parametrized with respect to the Gaussian, it is smaller at lorentz angle
double expNorm = CAMath::Sqrt(mDyExpNormA + mDyExpNormC * (snp - mLorentzAngle) * (snp - mLorentzAngle));
if (slope < convertAngleToDy(snp)) {
// left tail in this case, larger tail for large positive snp
likelihood += expNorm * CAMath::Exp((mDyExpA + mDyExpC * (snp - mLorentzAngle)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)));
}
else {
} else {
// right tail, larger tail for large negative snp
likelihood += expNorm * CAMath::Exp(-(mDyExpA - mDyExpC * (snp - mLorentzAngle)) * (slope - GPUTRDRecoParam::convertAngleToDy(snp)));
}

// Normalized such that likelihood is 1 when track angle and tracklet slope agree, and always lower than 1 otherwise (such that -log(likelihood) is always positive)
likelihood /= (1.f + expNorm);
return likelihood;
}

GPUdi() double GPUTRDRecoParam::getZLikelihood(float deltaZ, float padLength, float sigmaZtrk) const
GPUdi() double GPUTRDRecoParam::getZLikelihood(float deltaZ, float padLength, float sigmaZtrk) const
{
// logistic function as approximation of convolution between uniform tracklet Z and gaussian track, depends on sqrt(3)/pi times the track resolution
// normalized so that maximum is 1
double lmax = 1.f/(1.f + CAMath::Exp(- 0.5f * padLength / 0.5513f / sigmaZtrk)) - 1.f/(1.f + CAMath::Exp(0.5f * padLength / 0.5513f / sigmaZtrk));
double lmax = 1.f / (1.f + CAMath::Exp(-0.5f * padLength / 0.5513f / sigmaZtrk)) - 1.f / (1.f + CAMath::Exp(0.5f * padLength / 0.5513f / sigmaZtrk));
return 1.f / lmax / (1.f + CAMath::Exp((deltaZ - 0.5f * padLength) / 0.5513f / sigmaZtrk)) - 1.f / lmax / (1.f + CAMath::Exp((deltaZ + 0.5f * padLength) / 0.5513f / sigmaZtrk));
}

Expand Down
42 changes: 20 additions & 22 deletions GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -635,14 +635,13 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
if (probBC > maxProb) {
maxProb = probBC;
yCorrPileUp = -slopeFactor * deltaBC;
zShiftTrkPileUp = - deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
}
}
if (sumProb > 1e-6f) {
yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
}
}


// number of tracklets within the chamber is the current TRD occupancy estimator
int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet];
Expand All @@ -652,13 +651,13 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
float zPosCorr = spacePoints[trkltIdx].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl();
float yPosCorr = spacePoints[trkltIdx].getY() - tiltCorr + yCorrPileUp;
zPosCorr -= zShiftTrk + zShiftTrkPileUp; // shift tracklet instead of track in order to avoid having to do a re-fit for each collision

// Correction of y position based on angular pull
if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
yPosCorr += corrPull;
}

float deltaY = yPosCorr - projY;
float deltaZ = zPosCorr - projZ;

Expand All @@ -677,18 +676,18 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
trkltCovTmpWithDy[0] += trkWork->getSigmaY2();
trkltCovTmpWithDy[1] += trkWork->getSigmaZY();
trkltCovTmpWithDy[2] += trkWork->getSigmaZ2();

if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
// In this case the correlation between tracklet y and tracklet dy is already corrected for so we can cancel it
trkltCovTmpWithDy[3] = 0.;
trkltCovTmpWithDy[4] = 0.;
}

// We add the correlation between track y and track dy
trkltCovTmpWithDy[3] += trkWork->getSigmaSnpY() * mGeo->GetCdrHght();
trkltCovTmpWithDy[4] += trkWork->getSigmaSnpZ() * mGeo->GetCdrHght();
// keep these parameters before matrix inversion

// keep these parameters before matrix inversion
float sigmaZ2 = trkltCovTmpWithDy[2];
float sigmaDy2 = trkltCovTmpWithDy[5];

Expand All @@ -698,13 +697,15 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
if (Param().rec.trd.addDeflectionInChi2 == 2 || Param().rec.trd.addDeflectionInChi2 == 3) {
// In this case we take into account the full likelihood, so we replace (deltaDy/sigmaDy)^2 by -2*ln(likelihood), which is the same in the default Gaussian case
double likelihood = mRecoParam->getDyLikelihood(trkWork->getSnp(), spacePoints[trkltIdx].getDy() + dyTiltCorr, nTrackletsChamber);
if (likelihood < 1e-6f) continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation
if (likelihood < 1e-6f)
continue; // likelihood of 1e-6 is equivalent to 5 sigma deviation, so we can safely cut it to avoid numerical instability in log calculation
deltaDy = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaDy2) * (deltaDy > 0.f ? 1.f : -1.f);
}
if (Param().rec.trd.addDeflectionInChi2 == 3) {
// We do the same for deltaZ
double likelihood = mRecoParam->getZLikelihood(deltaZ, pad->GetRowSize(tracklets[trkltIdx].GetZbin()), CAMath::Sqrt(trkWork->getSigmaZ2()));
if (likelihood < 1e-6f) continue;
if (likelihood < 1e-6f)
continue;
deltaZ = CAMath::Sqrt(-2.f * CAMath::Log(likelihood) * sigmaZ2) * (deltaZ > 0.f ? 1.f : -1.f);
}
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;
Expand All @@ -724,8 +725,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
Hypothesis hypoNoUpdate(trkWork->getNlayersFindable(), iCandidate, -1, trkWork->getChi2() + Param().rec.trd.penaltyChi2);
InsertHypothesis(hypoNoUpdate, nCurrHypothesis, hypothesisIdxOffset);
isOK = true;
} // end candidate loop

} // end candidate loop

mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->getChi2(), iLayer); // only meaningful for ONE candidate!!!
mDebug->SetRoad(roadY, roadZ, iLayer); // only meaningful for ONE candidate
Expand Down Expand Up @@ -791,9 +791,9 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
pad = mGeo->GetPadPlane(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector());
float tiltCorrUp = tilt * (spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() - trkWork->getZ());
float dyTiltCorr = tilt * trkWork->getTgl() * mGeo->GetCdrHght();

float yPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getY() - tiltCorrUp;

float zPosCorrUp = spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getZ() + mRecoParam->getZCorrCoeffNRC() * trkWork->getTgl();
zPosCorrUp -= zShiftTrk;
float padLength = pad->GetRowSize(tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin());
Expand Down Expand Up @@ -825,21 +825,21 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
if (probBC > maxProb) {
maxProb = probBC;
yCorrPileUp = -slopeFactor * deltaBC;
zShiftTrkPileUp = - deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
zShiftTrkPileUp = -deltaBC * o2::constants::lhc::LHCBunchSpacingMUS * mTPCVdrift * mTrackAttribs[iTrk].mSide;
}
}
if (sumProb > 1e-6f) {
yAddErrPileUp2 = sumCorr2 / sumProb - 2 * yCorrPileUp * sumCorr / sumProb + yCorrPileUp * yCorrPileUp;
}
}

zPosCorrUp -= zShiftTrkPileUp;
yPosCorrUp += yCorrPileUp;

const auto currDet = tracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector();
int nTrackletsChamber = mTrackletIndexArray[trkltIdxOffset + currDet + 1] - mTrackletIndexArray[trkltIdxOffset + currDet];
float angularPull = GetAngularPull(spacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].getDy() + dyTiltCorr, trkWork->getSnp(), nTrackletsChamber);

// Correction of y position based on angular pull
if (Param().rec.trd.useAngularPull == 3 || Param().rec.trd.useAngularPull == 4) {
float corrPull = -angularPull * mRecoParam->getCorrYDy(trkWork->getSnp());
Expand Down Expand Up @@ -938,8 +938,6 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
return false;
}
} // end layer loop



// --------------------------------------------------------------------------------
// add some debug information (compare labels of attached tracklets to track label)
Expand Down Expand Up @@ -1288,7 +1286,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::IsGeoFindable(const TRDTRK* t, const
if (!mGeo->ChamberInGeometry(det)) {
return false;
}

// reject tracks in bad chambers
if (mChamberStatus[det]) {
return false;
Expand All @@ -1310,12 +1308,12 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::IsGeoFindable(const TRDTRK* t, const
if (!((zTrk > zMin + epsZ) && (zTrk < zMax - epsZ))) {
return false;
}

// reject tracks if the corresponding pad and neighboring pads (to take into account track uncertainty and charge sharing with neighbours) are masked
int32_t padrow = pp->GetPadRowNumber(zTrk);
int32_t padcol = pp->GetPadColNumber(t->getY());
int32_t idxPad = det * kNPadColumns * kNPadRows + padcol * kNPadRows + padrow;
if (padrow != -1 && padcol != -1 && mPadStatus[idxPad] && (padrow == 0 || mPadStatus[idxPad - kNPadRows]) && (padrow == kNPadRows-1 || mPadStatus[idxPad + kNPadRows])) {
if (padrow != -1 && padcol != -1 && mPadStatus[idxPad] && (padrow == 0 || mPadStatus[idxPad - kNPadRows]) && (padrow == kNPadRows - 1 || mPadStatus[idxPad + kNPadRows])) {
return false;
}

Expand Down
2 changes: 1 addition & 1 deletion GPU/GPUTracking/TRDTracking/GPUTRDTracker.h
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ class GPUTRDTracker_t : public GPUProcessor
kNChambers = 540,
kNPadRows = 16,
kNPadColumns = 144,
kNPads = 1244160};
kNPads = 1244160 };

struct HelperTrackAttributes {
// additional TRD track attributes which are transient
Expand Down
Loading